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Preface 


The  purpose  of  this  study  is  to  develop  a  code  based  on 
the  Approximate  Navier  Stokes  (ANS)  equations,  to  allow  flow 
computation  at  low  Reynolds  Number  (less  then  500,000)  about 
highly  cambered  airfoils.  The  camber  may  either  be  the 
airfoil  characteristic  or  created  by  deflecting  portions  of 
the  Leading  edge  (Slat)  or  Trailing  edge  (Flap),  or  both. 

This  study  specifically  looks  at  the  highly  cambered 
Wortman  FX  63-137  airfoil,  fitted  with  flaps  and  slats  at 
various  combinations  of  deflections  and  angles  of  attack. 

The  flow  has  been  computed  using  the  ANS  equations  in  delta 
form  of  the  Vorticity  Stream  Function  equations.  Qood 
agreement  with  experimental  results  has  been  obtained  for 
small  angles  of  attack. 

The  significance  of  this  study  is  based  on  the  fact  that 
recently  there  has  been  an  increased  application  of  low 
Reynolds  Number  aerodynamics,  Remotely  Piloted  Vehicles 
(RPV's),  Helicopter  rotors  and  turbine  blades  are  only  few  of 
such  applications.  This,  combined  with  the  dearth  of 
knowledge  in  this  area  has  created  an  environment  where 
detailed  information  of  the  flow's  structure  is  needed.  It 
is  hoped  that  this  study  will  help  fill  a  portion  of  this 
knowledge  gap  by  providing  to  the  practicing  aerodynamicist 
and  experimenters,  the  tool  to  allow  numerical  flow 
prediction,  and  a  reference  data  set  on  the  wortman  airfoil, 
which  is  contained  herein. 
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Abstract 


/ 

The  purpose  of  this  study  is  to  develop  a  code  based  on 
the  Approximate  Navier  Stokes  (ANS)  equations  (in  the 
Vorticity  Stream  Function  delta  form).  The  Wortman  FX  63-137 
airfoil  fitted  with  leading  and  trailing  edge  devices  has 
been  analyzed  at  low  Reynolds  Number  (100,000)  at  various 
angles  of  attack  and  various  deflection  angles,  results  are 
compared  to  the  experimental  data.  The  agreement  is  very 
good  at  small  angles  of  attack.  However,  at  large  angles  of 
attack  the  disparity  is  larger,  and  improvements  have  been 
recommended  which  will  resolve  this  disparity.  Overall  the 
present  scheme  produces  very  reasonable  results,  with  good 
repeatability  and  fast  convergence,  and  has  the  potential  of 
being  developed  into  an  effective  design  tool.  Tkei^':  .  /V  rl 
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NUMERICAL  STUDY  OF  THE  INFLUENCE  OF  LEADING  AND  TRAILING 


EDGE  FLAPS  ON  THE  PERFORMANCE  OF  AIRFOIL 


I.  Introduction 


The  aim  of  every  aerodynamicist  is  to  be  able  to  predict 
the  aerodynamic  character ist ics  of  a  given  design.  Accurate 
prediction  of  these  parameters  is  very  important  as  they  form 
the  basis  for  the  designs  of  the  aircraft's  control  and 
propulsion  systems,  and  hence  effect  the  safety  of  the 
vehicle. 

The  tools  available  to  an  aerodynamicist  are  the  wind 
tunnels  and  computational  methods.  To  date,  most  aircraft 
designs  depend  on  wind  tunnel  testing  and  empirical  data, 
but,  lately,  some  designs  have  also  been  developed  primarily 
through  the  use  of  computational  methods.  This  has  largely 
been  possible  due  to  the  enormous  growth  of  research  and 
success  in  the  area  of  Computational  Fluid  Dynamics  (CFD) . 
Computational  methods  are  more  economical  and  faster  compared 
to  wind  tunnel  testing  and  give  much  more  detailed 
information  on  flow  characteristics.  Geometric  changes  and 
their  effects  on  flow  can  also  be  seen  very  quickly, 
resulting  in  reliable  design  decisions. 

Computational  methods  can  be  broken  down  into  two 
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distinct  techniques,  finite  element  and  finite  difference 
algorithms.  This  research  is  based  on  the  finite  difference 
approach  and  is  an  attempt  to  develop  a  fast  and  efficient 
algorithm  for  computing  flow  about  arbitrary  bodies 
(primarily  airfoils),  using  the  approximate  Navier-Stokes 
(ANS)  equations.  These  equations  are  not  as  accurate  as  the 
Navier-Stokes  (NS),  but  the  advantages  are  that  the  ANS 
equations  offer  significantly  reduced  computational  effort 
and  storage  compared  to  the  complete  Navier  Stokes  equations. 

Purpose  and  Scope 

The  basic  purpose  of  this  research  is  to  develop  a  code 
which  would  allow  flow  computation  at  low  Reynolds  Number 
(less  then  500,000)  about  highly  cambered  airfoil  sections. 
The  high  camber  is  usually  produced  by  installing  leading 
edge  (Slats)  and  trailing  edge  (Flaps)  devices  to  the  basic 
section.  These  devices  can  then  be  rotated  to  effectively 
change  the  shape  (camber)  and  hence  the  aerodynamic 
characteristics  of  the  airfoil. 

Significance.  Motivation  for  this  research  is  the  fact 
that  extensive  experimental  and  empirical  data  is  available 
for  airfoils  with  flaps  and  slats  at  high  Reynolds  number, 
but  very  little  is  known  about  the  aerodynamic  behavior  of 
highly  cambered  airfoils  in  low  Reynolds  number  f!cw.  The 
second  factor  was  the  ready  availability  of  experimental  data 
for  the  Wortman  FX  63-137  in  various  configurations  and  at 
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varying  Reynolds  number  from  the  University  of  Notre  Dame 
[16]  and  the  College  of  Aeronautics,  Cranfield,  England  [28]. 
The  third  and  most  important  factor,  is  that  recently  there 
has  been  recently  increased  application  of  low  Reynold  number 
aerodynamics,  remotely  piloted  vehicles  (RPV),  helicopter 
rotors,  turbine  blades  and  wind  turbines  are  only  few  of 
such  applications.  These  are  all  character ized  as  having 
either  short  chord  length  or  low  operating  speed  or  both.  A 
diagram  showing  the  range  of  Reynolds  Numbers  for  several 
applications  appears  below  as  Figure  1.  Therefore,  it  is 
very  natural  to  anticipate  a  vast  need  for  performance 
character istics  in  the  low  Reynolds  number  regime,  and  how 
these  may  be  enhanced  using  various  techniques. 


Before  the  issue  of  performance  enhancement  can  be 
discussed,  it  would  be  appropriate  to  highlight  a  typical 
character ist ic  of  airfoils  in  low  Reynolds  Number  flow,  the 
formation  of  separation  bubbles  and  their  role  in  determining 
the  flow  pattern  and  the  performance  characteristics  of  the 
airfoi 1 . 

Extensive  research  has  been  carried  out  on  separation 
bubbles,  and  it  has  been  determined  that  the  significant 
parameters  affecting  the  formation  of  separation  bubbles  are 
the  angle  of  attack,  profile  geometry,  Reynolds  Number  and 
free  stream  turbulence  [16:3].  The  separation  bubbles  are 
usually  classified  as  either  long  or  short,  Arena  and  Muller 
[3]  summarize  that  long  bubbles  tend  to  decrease  lift  due  to 
lower  suction  peak,  while  short  bubbles  have  negligible 
effect  on  lift  or  suction  peak. 

The  separation  bubble  forms  as  the  fluid  moving  over  the 
airfoil  encounters  an  adverse  pressure  gradient.  If  the  flow 
does  not  possess  sufficient  momentum  to  overcome  this 
gradient,  the  laminar  boundary  layer  separates  from  the 
surface  and  then  undergoes  transition  to  turbulent  flow.  It 
then  gains  momentum  by  entraining  fluid  from  the  free  stream, 
and  if  the  fluid  gains  sufficient  momentum  to  overcome  the 
gradient,  It  will  reattach  to  the  airfoil  surface.  Hence 
the  process  of  separation  and  reattachment,  creates  a 
reverse  flow  region  over  the  airfoil  surface,  usually  known 
as  the  separation  bubble.  On  the  other  hand,  if  the  flow  is 
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unable  to  gain  sufficient  momentum  to  reattach,  then  it 
remains  separated.  This  type  of  separation  is  known  to  yield 
an  open  wake  of  an  airfoil  and  is  the  primary  cause  of  stall. 

The  performance  of  an  airfoil  may  be  enhanced  by 
management  of  the  boundary  layer  using  any  one,  or 
combination,  of  the  following  techniques. 

(a)  Attachment  of  leading  and  or  trailing  edge  devices 

to  the  airfoi 1 . 

(b)  Use  of  surface  roughness. 

(c)  Suction  and  blowing. 

Because  this  study  is  concerned  primarily  with  the 
effects  of  leading  and  trailing  edge  devices,  it  would  be 
appropriate  to  draw  on  general  conclusions  based  on  previous 
work  done  in  this  area. 

Leading  Edge  Device  (Slat).  The  primary  purpose  of 
these  devices  is  to  avoid  or  delay  leading  edge  separation, 
particularly  at  low  Reynolds  number.  Positive  deflections  of 
the  leading  edge  (downwards)  causes  an  effective  reduction 
of  the  angle  of  attack  relative  to  the  airfoil's  leading 
edge,  and  leads  to  downstream  shift  of  the  separation  point. 

Trailing  Edge  Device  (Flap).  Positive  deflection  of  the 
flap  (downward)  causes  an  increase  in  pressure  acting  on  the 
lower  surface  and  a  decrease  of  pressure  on  the  upper 
surface.  This  produces  an  increase  in  the  circulation  around 
the  airfoil,  thereby  increasing  lift,  but  only  at  the 
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penalty  of  increased  drag.  The  increase  in  drag  is  due  to  an 
increase  in  the  effective  angle  of  attack  of  the  trailing 
section  of  the  airfoil,  and  low  pressure  in  the  wake 
downstream  of  the  flap.  The  overall  effect  of  the  flap  is  to 
shift  the  lift  curve,  without  changing  the  slope.  This  may 
be  predicted  by  the  following  relation,  derived  using  the 
thin  airfoil  theory  [12:492-517,14]. 

C^=  2n(a+2z) 

where  a  is  the  angle  of  attack  and  z  is  the  maximum  camber 

It  can  be  seen  from  this  relation  that  the  lift 
coefficient  will  increase  linearly  as  the  camber  is 
increased.  It  should  be  mentioned  that  this  relation  is  only 
applicable  to  thin  airfoils  in  ideal  flow,  but  is  not 
applicable  to  highly  cambered  airfoils  in  viscous  flow. 

Viscous  effects  are  generally  nonlinear  and  therefore 
cannot  be  approximated  by  such  simple  linear  relations. 

Hence,  to  be  able  to  predict  complex  viscous  flow  patterns 
effectively,  we  have  to  turn  to  mathematical  models  using 
either  the  Navier-Stokes  or  the  approximate  Navier-Stokes 
(ANS)  equations.  The  main  thrust  of  this  study  is  to  use  the 
ANS  model  to  approximate  the  flow  pattern  and  the  performance 
characteristics  of  the  Wortman  FX  63-137  airfoil,  which  has 
been  appropriately  modified  with  slats  extending  to  13%  and 
flaps  measuring  25%  chord  (a  Fortran  program  named  "AIRFOIL" 
has  been  developed  to  carry  out  deflections  of  the 
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airfoil).  To  be  able  to  analyze  this  problem  the  code 
developed  by  Halim  [11]  has  been  extended  from  the  purely 
steady  state  laminar  symmetric  flow  solver  to  the  present 
version.  The  present  version  is  flexible,  easy  to  use,  and 
can  handle  nonsymmetric  turbulent  flows.  Using  this  code  in 
either  steady  state  or  time  accurate  unsteady  solutions  can 
be  obtained  and  either  Neumann  or  Dirichlet  Boundary 
conditions  can  be  enforced.  Routines  have  also  been  added  to 
compute  pressure  distribution  and  the  global  parameters  lift 
(C^),  drag  (C^)  and  moment  (C^)  coefficients.  Also  as  a  part 
of  this  effort  an  elliptic  grid  generator  has  been  developed, 
which  allows  the  user  complete  control  over  all  aspects  of 
the  grid  generation  process.  A  special  wake  model  handles 
the  interaction  of  the  trailing  edge  deflection  and  the  angle 
of  attack  to  produce  a  C-type  grid  which  honors  these 
conditions . 

The  code  developed  in  this  study  is  based  on  the 
approximate  Navier-Stokes  equations  in  the  vorticity,  stream 
function  delta  form  and  has  been  validated  by  computing  and 
comparing  three  different  types  of  flows:  symmetric,  low 
Reynolds  Number  (Re=12,500),  laminar  flow  around  NACA-12 
airfoil,  turbulent  flow  on  a  finite  flat  plate  (Re=6. 5x10s) 
and  turbulent  flow  around  the  nonsymmetric  Wortman  FX  63-137 
airfoil  (Re=200 ,000) .  Each  of  these  test  cases  probes  one 
or  more  specific  areas  of  the  code,  which  allows  a 
systematic  validation  process.  The  results  obtained  in  each 
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of  the  three  cases  are  in  very  good  agreement  with  previous 
data. 

Having  completed  the  validation  process,  flow  around 
the  Wortman  airfoil  was  computed  at  Re=100,000.  Because  of 
time  constraint,  one  flap  deflection  of  20  degrees,  one  slat 
deflection  of  5  degrees  was  analyzed  at  different  angles  of 
attack,  in  addition  to  the  case  of  no  deflection.  Results 
obtained  are  compared  to  previous  data  [16,26],  agreement 
with  experimental  data  is  good  at  lower  angles  of  attack, 
with  under  prediction  of  C-  and  over  prediction  of  C.  and  C 

*C  CL  ^ 

at  higher  angles  of  attack  (-4  <  a  >  8  deg). 

Based  on  the  present  results,  an  attempt  has  been  made 
to  look  at  the  flow  structure  more  closely,  which  is  not 
alwayspossible  experimentally  (the  measurement  process  may 
alter  the  flow  pattern  and  gives  rise  to  un-certainty  [7]). 

It  is  hoped  that  the  insight  so  gained  will  benefit  future 
experiments  and  vehicle  designs. 

The  present  system  of  predicting  flow  can  be  improved  so 
as  to  allow  good  prediction  capabilities  by  using  different 
turbulence  models  and  some  minor  improvements  in  the 
algorithm.  Overall,  it  can  be  said  that  the  present  code 
has  good  potential  of  developing  into  an  effective  design 
tool  and  this  study  is  only  the  first  step  in  this  direction. 

The  following  chapters  discuss  at  length,  details  of 
the  mathematical  model,  the  solution  algorithm,  the 
validation  process  and  the  results  for  the  Wortman  airfoil, 
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finally  suggestions  have  been  made  to  enhance  the  present 
work  and  specific  areas  identified  in  which  an  in-depth 
analysis  is  required. 
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I  I  Analysis 


In  this  chapter,  the  mathematical  model  used  in  this 
study  will  be  developed,  starting  from  the  fundamental 
equations  for  viscous  flow  and  building  up  to  the  model 
finally  used  in  finite  difference  form. 


Basic  Equations 

The  fundamental  equations  governing  two-dimensional 
incompressible  flow  of  a  fluid  with  no  body  forces  and 
constant  properties  are  the  Navier  Stokes  (NS)  equations 
[20:65]. 


continuity 


du  dv  _ 
3x  dy  " 


0 


x-mom. 


du 

dx 


du  _  2±_  dP 
dy  p  5x 


+ 


O) 


y-mom. 


dv 

dx 


dv  _  -1_  dP 
vdy  p  ay 


It  should  be  noted  that  equations  (1)  are  in  dimensional 
form  and  viscosity  has  not  been  assumed  constant.  This  has 
been  done  to  facilitate  simulation  of  turbulent  flow,  using 
the  concept  of  eddy  viscosity. 


Non-dimensional ization 

The  above  equations  (1)  can  be  non-dimensional ized  using 
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the  following  normalization  definitions. 


x 


P  = 


p*u*  2 
K®  ® 


p  = 


A 

)_ 

A 


Re  = 


A  A 

p  U  C 
K®  ® 


® 


(2) 


where  c  is  the  characteristic  length  (airfoil  chord:  c  =  1) 
and  subscript  ®  indicates  free  stream  properties. 


In  the  above  relations  starred  (♦)  variables  represent 
dimensional  quantities  and  non-starred  variables  represent 
the  nondimensional  variables.  Equations  (2)  take  the 
following  form  after  normalization. 


du  dv  _ 
<5x  dy 


0 


(3b) 


du  *  .&L  4.  xjdu  -  dP  4.  1  9  f  9ul  ,  3  f  diA  » 

sFE  +  u5x  +  vWy  ~  !Fx  +  ~Re  +  (3a) 

dv  4.  4  „9v  -  9P  .  1  a  [  dv'l  .  d  C  dv|  x 

5t+ua5+v5y-^y  +  “R5  a?[^J  +  Wy[^y)  (3c) 


Stream  function  and  Vorticity  transport  equations 

The  above  non-dimensional ized  equations  (3)  can  be  cast 
into  the  stream  function  and  vorticity  transport  equation  by 
differentiating  equation  (3b)  with  respect  to  y  and 
differentiating  equation  (3c)  with  respect  to  x.  Subtracting 
the  two  resulting  equations  from  each  other  eliminates  the 
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pressure  term,  and  the  rest  of  the  terms  can  be  simplified 
using  the  following  definitions  for  vorticity  and  stream 
function,  (for  more  details,  see  Appendix  A) 

(vorticity)  o  =  §S  +  fp  (4a) 

(stream  function)  u  =  ^  &  v  =  -  (4b) 

where  is  the  stream  function  and  <■)  the  vortic 

The  Naviei — Stokes  equations  in  vorticity  and  stream 
function  form  are 

(vorticity  equation) 

“t  +  V*  -  *xuy  =  H  XX  ♦  H  yy 

2#ixy*xy  +  ^xxvyy  +  *VyVxx  (5a) 

(stream  function  equation) 

v  +  v  =  (<)  (5b) 

rxx  wyy  v  7 

here  subscripts  indicate  differentiation 

In  the  above  set  of  equations  (5),  the  NS  equations  in 
primitive  variables,  (3b)  &  (3c),  transform  to  the  vorticity 
transport  equation  (5a).  The  continuity  equation  (3a)  is 
satisfied  identically  by  the  definition  of  stream  function 
equation  (3a).  Thus  equations  (5)  represent  two  unknown 
variables  (■>  and  ?  in  two  equations,  forming  a  closed  set, 
provided  p  is  known  from  the  turbulence  modeling.  The 
solution  of  (5)  can  be  obtained  by  solving  them 
simultaneously,  or  through  a  procedure,  that  will  be 
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discussed  in  chapter  4  in  detail. 


Coordinate  transf ormat ion 

The  stream  function  and  vorticity  transport  equations 
are  in  the  Cartesian  coordinates  and  therefore  can  be  applied 
easily  to  only  rectangular  domains.  To  facilitate 
application  to  irregular  boundaries,  a  transformat  ion  is 
required  that  would  map  the  equations  and  the  boundary 
conditions  from  the  irregular,  physical  domain  to  a  body 
conformal  coordinate  system. 

The  physical  domain  represents  the  actual  geometry  of 
the  problem  and  may  be  of  any  shape,  as  shown  in  Figure  2. 
Let  x  and  y  be  the  co-ordinate  variables  in  the  physical 
domain  and  rj»  €  be  the  co-ordinate  variables  in  the 
computational  domain  as  shown  in  Figure  2.  Then  if  the 
following  transformation  exists; 

5  =  C(x,y) 
n  =  n (x,y) 

The  vorticity  stream  function  equations  can  be 
transformed  into  the  computational  domain  using  the  chain 
rule  of  partial  differentiation.  This  results  in  the 
co-ordinate  transformation  relations  (25:122-133  )  given  in 
Appendix  B. 

Parabol ization 

The  complete  Navier  Stokes  equations  are  the  most 
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accurate  for  any  flow  situation,  but  these  are  very 
difficult  to  solve  in  their  complete  form.  In  general,  a 
very  large  amount  of  computer  storage  and  computational  time 
is  necessary  to  obtain  a  solution  of  these  equations  [2:420]. 
However  for  some  situations,  assumptions  may  be  made  that 
will  reduce  these  equations  to  a  simpler  form.  Boundary- layer 
(BL)  equations  and  the  Approximate  Navier  Stokes  (ANS) 
equations  fall  into  that  category.  Following  are  some  of  the 
advantages  of  using  the  ANS  equations  [2:420]. 

(a)  There  are  fewer  terms  in  the  equations,  which  leads 
to  reduction  in  computation. 

(b)  The  equations  are  parabolic  in  the  streamwise 
direction,  so  that  their  solution  can  be  obtained  using 
spatial-marching  techniques.  This  result's  in 
substantial  savings  of  computational  time  and  storage 
requirements. 

(c)  Unlike  the  boundary- layer  equations,  which  are 
singular  at  the  point  of  separation,  the  ANS  equations 
can  easily  handle  reverse  flow. 

The  derivation  of  the  ANS  equations  from  the  NS 
equations  is  not  as  rigorous  as  the  derivation  for 
Boundary- 1 ayer  equations.  However,  the  following  assumptions 
are  usually  made  to  arrive  at  the  ANS  form  of  the  NS 
equations . 

(a)  The  normal  pressure  gradient  term  is  retained. 
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(b)  The  second  derivatives  with  respect  to  the 
streamwise  direction  are  omitted. 


The  assumptions  given  above  have  also  been  applied  in 
this  work,  and  equations  (4)  reduce  to  the  following  set 
(after  the  transf ormat ion  into  general  coordinate  system  and 
application  of  ANS  assumptions). 

stream  function 

ff(~X  H  "  “j"  In)  +in(“J"  §7  “  “5“  H)  +Jw=0  (6a) 

vorticity  transport 

J H ♦  S5 5! - 3 -Ji  +  <«*> 
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Ill  Numerical  solution 


In  this  chapter,  numerical  grid  generation  is 
discussed,  followed  by  the  numerical  method  for  solving  the 
vorticity  stream  function  equations,  the  differencing 
scheme,  the  boundary  conditions,  the  turbulence  closure. 

The  transition  from  laminar  to  turbulent  flow,  the 
convergence  criterion,  computation  of  pressure  and  force 
coefficients. 

Grid  Generation 

As  mentioned  earlier  the  purpose  of  numerical  grid 
generation  is  to  remove  the  problem  of  boundary  shape  from 
finite-difference  methods.  Numerically  generated  grids  allow 
implementation  of  numerical  algorithms  on  a  rectangular 
computational  domain,  regardless  of  the  shape  and 
configuration  of  the  physical  region  [25:v]. 

Numerical  grid  generation  is  not  in  the  scope  of  this 
thesis.  But  efforts  were  made  to  develop  a  grid  generator, 
primarily  because  a  grid  was  required  to  solve  the  present 
problem.  Secondly  no  available  grid  generator  could  handle 
the  specifics  of  this  study,  these  being. 

(a)  The  ability  to  numerically  regenerate  the  surface 
of  the  airfoil  so  as  to  simulate  leading  and  trailing 
edge  deflections. 

(b)  The  grid  generator  should  also  be  able  to  handle 
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the  combined  effect  of  trailing  edge  deflections  and 
angle  of  attack  on  the  wake  branch  cut  (C-grid). 

(c)  Capability  to  control  spacing  and  clustering  of 
points  anywhere  in  the  physical  field  and  on  the 
boundaries . 

Airfoil  Deflections.  The  solution  of  problem  (a) 
involves  numerical  rotation  ( transformat ion)  of  airfoil 
coordinates  about  rotation  centers  based  upon  given 
deflection  data.  A  code  named  "AIRFOIL"  was  developed 
(Fortran  listing  attached  as  Appendix  C),  that  takes  as  input 
the  airfoil  coordinates,  lengths/deflections  of  the  leading 
and  trailing  edge  devices  and  regenerates  the  surface 
satisfying  the  given  conditions,  Figure  3  represents  a 
comparison  of  the  basic  Wortman  section  before  and  after 
rotating  the  leading  edge  by  15°  and  trailing  edge  by  20° 
using  the  program  AIRFOIL. 

The  code  Airfoil  also  has  the  ability  to  smooth 
coordinate  data,  so  that  C1  (first  derivative)  continuity  of 
the  metrics  is  ensured  on  the  surface.  The  smoothing  of 
coordinates  is  implemented  with  an  IMSL  routine  that  uses 
cubic  splines  and  least  squares  to  minimize  metric  variation. 
The  addition  of  smoothing  was  necessitated  for  to  the 
following  reasons. 

(a)  After  rotation  of  the  airfoil  coordinates  some  of 
the  points  have  to  be  removed  or  altered  due  to  overlap  or 
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(a)  No  deflections 


(b)  deflections:  Leading  edge  15°  trailing  edge  20° 

Figure  3  Comparison  of  un-defl acted  and  deflected  Wortman 
FX63-137  Airfoil 

separation.  A  first  cut  to  this  approximation  is  done  as 
part  of  the  rotation  itself,  but  this  does  not  ensure  C1 
continuity,  therefore  a  need  arises  to  apply  mathematical 
smoothing. 

(b)  The  coordinate  of  the  Wortman  FX  S3-137  airfoil 
obtained  from  references  [16:110]  and  [5:688],  are  not  C1 
smooth  which  caused  problems  with  solution  convergence  and 
inaccuracy,  and  forced  unwanted  spikes  and  noise  into  the 
computed  data.  Thus,  smoothing  of  the  airfoil  coordinates 
was  necessary. 

Figure  4  shows  that  smoothing  has  considerably  improved 
the  C1  variation  along  the  surface  (though  it  could 
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have  been  better),  without  significantly  changing  the 
geometry  of  the  airfoil.  Figure  5,  Appendix  D  tabulates  the 
coordinates  before  and  after  smoothing.. 

Wake  Prof i 1 e .  The  problem  of  approximating  the  wake 
profile  or  the  branch  cut  for  the  C-grid  has  been  solved  by 
piecing  together  three  polynomials.  The  first  segment  is  a 
sixth-order  polynomial,  which  approximates  the  start  of  the 
profile  based  on  mean  trailing  edge  deflection  and  location. 
The  first  patch  terminates  0.6c  downstream  of  the  trailing 
edge,  where  the  final  wake  angle  is  enforced.  The  wake 
angle  can  be  either  a  user  input  or  the  angle  of  attack 
(default) . 

Figure  6  shows  a  sketch  of  the  wake  profile. 

Coefficients  a  through  a  of  the  polynomial  approximating 
0  6 

this  curve  are  found  by  applying  the  following  conditions, 
(a)  x  =  0  y  =  0 
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The  second  patch  is  nearly  the  same  as  the  sixth  order 
polynomial  described  above.  The  only  differences  being  in 
the  values  of  x^  and  x ,  which  are  so  chosen  as  to  be 
optimum  for  this  patch,  and  would  allow  y  variation  without 


oscillations  (characteristic  of  higher-order  polynomials). 

This  patch  picks  the  slope  error  at  the  end  of 
first  patch  and  blends  it  down  to  zero.  Thus  by  the  end  of 
this  patch,  the  profile  has  acquired  the  required  slope. 

The  third  and  final  patch  linearly  extrapolates  the 
co-ordinates  till  the  end  of  the  branch  cut  based  on  the 
slope  information  provided  to  it  by  the  second  patch. 

The  overall  problem  of  grid  generation  (c)  poses  the 
most  stringent  requirement,  and,  in  the  opinion  of  the 
author,  is  best  solved  by  an  elliptic  grid  generation  scheme, 
such  a  scheme  has  the  inherent  property  of  generating  very 
smooth  grids  and  provides  complete  control  over  grid  point 
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distribution.  A  core  program  was  thus  acquired  from  Lt. 
Amdhal  (AFWAL )  [1],  that  solves  equations  (7)  to  generate  a 
grid  (based  on  the  boundary  point  spacing). 

oxEE  ’  2e*5r,  *  >xnn  =  '-'2<'>x5  ♦  %> 

oy«  '  JSy?n  +  iyrm  =  'j2(Pyt  *  %’  (7) 

where  a,  0,  y ,  J  are  defined  in  Appendix  B,  and  P  and  Q  are 
terms  used  to  control  point  spacing  in  the  interior  of  the 
domain. 

Because  the  core  program  required  boundary  definition 
data  to  generate  grids,  a  program  named  "BOUNDARY"  was 
developed.  This  program  is  interactive  and  allows  the  user 
complete  control  over  point  distribution  and  spacing  on  all 
four  boundaries  of  the  C-grid  see  Figure  7. 

This  program  borrows  some  stretching  and  curve  fitting 
routines  from  a  hyperbol ic-gr id  generator  developed  by  Kinsey 
and  Barth,  AFWAL  [13],  and  employs  the  Vinokur  stretching 
function  [26]  for  point  distribution  and  spacing  variation. 
The  hyperbolic-grid  generator  is  based  upon  the  algorithm 
developed  by  Steger  and  Chaussee  [22]. 

The  elliptic-grid  generator  thus  creates  grids  based  on 
boundary  point  and  spacing  output  of  the  program  Boundary. 

An  example  of  a  grid  so  generated,  is  given  as  Figure  8. 

In  practice  however,  it  was  found  that  elliptic  grid 
generation  does  indeed  provide  complete  control  over  all 
aspects  of  the  grid,  but  has  the  serious  short-coming  of 
being  computationally  expensive  .  Also,  convergence  is  not 
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Figure  8  C-grid  created  by  the  Elliptic  grid  generator 


always  ensured  and  if  it  does  occur  it  can  take  a  very  long 
time  to  converge  to  an  acceptable  tolerance.  Therefore,  the 
grids  used  in  this  study  were  generated  using  a  modified 
version  of  the  hyperbolic  grid  generator,  described  earlier 
in  this  section.  The  modifications  were  primarily  restricted 
to  generation  of  the  wake  profile. 


1  n 


X  c 

Figure  9  C-grid  generated  using  Hyperbolic  grid  generator 

The  advantages  of  this  grid-generation  system  are  that 
it  ia  unconditionally  stable  and  produces  smooth  and  nearly 
orthogonal  grids.  Furthermore  it  is  very  inexpensive  to  use 
and  generates  grids  interactively  in  a  very  short  time.  The 
disadvantage  is  that  very  little  control  is  available  on 
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spacing  and  distributions  of  grid  lines  in  the  marching 
direction  (normal  to  the  inner  boundary).  The  theory  and 
program  details  of  the  hyperbol ic-gr id  generation  system  are 
discussed  at  length  in  Ref.  [13].  Figure  9  is  an  example  of 
a  grid  generated  hyperbol ical ly . 

Numerical  method 

The  vort icity-transport  equation,  (6a)  consists  of  an 

unsteady  term  w.  ,  the  advective  terms  w.u  and  w  «„  ,  and 

t  c  n  n  c 

the  viscous  diffusion  term  (uu)  .  This  equation  is  nonlinear 

m 

in  the  advective  terms,  and  must  be  linearized  before  any 
numerical  algorithm  can  be  applied  to  it.  Furthermore,  the 
equation  is  parabolic  in  time  and  in  the  streamwise  direction 
(5),  but  due  to  the  nonlinearity  must  be  solved  iteratively. 
On  the  other  hand,  the  the  stream  function  equation  is  linear 
and  elliptic  and  may  be  solved  using  either  direct  methods  or 
iterative  techniques. 

Basically  there  are  two  approaches  to  solve  the 
vort icity-stream  function  equations,  one  is  the  time 
dependent  approach  (unsteady  problem)  and  the  other  is  the 
steady-flow  approach  (steady  state  problem)  (18:107).  In 
this  study,  both  of  these  approaches  have  been  incorporated 
into  the  code  so  that  the  user  may  benefit  from  the 
advantages  of  either  scheme.  A  comparison  of  the  two 
approaches  follows. 

Time-dependent  approach.  The  Poisson  equation  for 


stream  function  is  iterated  to  convergence  at  each  time  step 
of  the  vort icity-transport  equation.  The  solution  is  marched 
in  time,  therefore  is  time  accurate.  This  approach  is  more 
stable  numerically  than  the  steady  state  approach,  but  is 
computationally  intensive. 

Steady-state  approach.  In  this  approach,  the  stream 
function  equation  is  not  iterated  to  convergence  at  each  time 
step.  Instead  both  equations  are  iterated  simultaneously, 
and,  if  convergence  takes  place,  steady  state  can  be  reached 
more  quickly  than  with  the  time-dependent  approach.  However, 
because  the  stream  function  equation  is  not  converged  at  each 
vorticity  iteration,  a  poor  initial  condition  can  lead  to 
non-linear  instability  (18:107-108). 

The  pseudo- 1  inear ired  vorticity  equation  is  given  below 
along  with  the  stream  function  equation.  These  equations  are 
expressed  for  a  general  curvilinear  coordinate  system  and  are 
in  delta  (6)  form  (see  Appendix  C  for  derivations). 

J  £e+  VV  V“e  ’  “n{»{  ‘  W 

=  -  ♦  -*;(-$-{*  *>)„)„  {8a) 

Hr*t  ‘  Hr  *„){+  (-5-  -  Hr  «*t)I1+  ■*“ 

=  '  ("X  "  Hr  *n)  '  (-5-  *n  •  Hr  *{)"  J“  (8b) 
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Discretization 


In  order  to  be  able  to  handle  arbitrary  geometries,  a 
second-order  accurate  (in  space)  conservative  differencing 
scheme  is  generated  for  the  stream  function  equation  (8b)  by 
integrating  the  equation  around  a  differential  element  in  the 
physical  domain,  see  Figure  (10). 


Figure  10  Differential  element  [11:3] 

The  resulting  system  of  algebraic  equations  is  solved 
iteratively  using  SLOR  (Successive  Line  Over  Relaxation). 
For  each  £-line  (constant-i)  in  the  computational  domain, 
the  finite-difference  equations  at  each  node  can  be  written 
as; [11: 3]. 
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A*6y  ♦  B*6|i  +  C  8*>* 

i  r>  j-i  J  wi, J  J  m. J*i 

+  D"&>  ♦  E'&j  +  F“6u  =  R* 

J  i»J-i  J  i,J  J  * ,  J ♦  1  J 


(9a) 


A°6f  +  B°6v  +  C°$¥> 
J  i.J-i  J  i.J  J 


+  D°&) 


+  E°6u  +  F°6u 


j  i.J-i  J  t.J  J  i.J+i 


=  R° 


(9b) 


where  superscript  denotes  a  coefficient  of  the 
stream  function  equation;  "©"  identifies  coefficient  of  the 
vorticity  equation;  ”j"  indicates  the  location  of  the  node  in 
the  redirection,  and  " 6 is  the  change  between  two 
successive  iterations  or  time  levels,  given  by 


vn*1=vn+6^n 


(10) 


Thus  along  each  (-line  a  tridiagonal  system  of  equations 

is  generated  from  the  stream  function  equation,  which  can  be 

solved  using  the  Thomas  algorithm  until  convergence  is 

reached  at  each  step  of  the  vorticity  equation.  Note  that 

the  coefficients  D*  ,  E  *  and  F*  would  be  zero  for  the  stream 

J  J  i 

function  equation  in  time-accurate  mode). 

For  the  vorticity  equation,  (8a),  central-difference 
approximations  are  used  for  all  terms  except  for  the  term 
(y>^5(i>^),  which  is  treated  as  an  upwind  difference.  In  the 
limit  of  the  steady  state  a  second  order  accurate  solution  is 
obtained  ref .  [11). 

Thus,  at  each  node  along  an  (-line,  the  correction 
equations  from  Newton  linearization  have  the  general  form 
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given  by  (9).  Equations  (9)  is  a  set  of,  a  block  (2x2) 
tridiagonal  system  of  equations,  and  is  solved  using  the 
modified  Thomas  algorithm  (simultaneous  solution  of  (9)  in 
steady  state  mode). 

Boundary  conditions 

Since  the  governing  equations  remains  the  same,  what 
distinguishes  one  flow  field  from  the  other  are  the  boundary 
conditions,  the  initial  conditions  and  of  course  flow 
parameters,  such  as  Reynolds  Number  (Re). 

The  specification  of  boundary  conditions  not  only 
affects  the  numerical  stability  of  a  computational  scheme, 
but  also  greatly  affects  the  numerical  accuracy  of  the 
finite-difference  solution.  The  boundary  conditions 
therefore  take  on  added  significance  and  must  be  based  upon 
physical  phenomenon. 

In  this  study,  boundary  conditions  are  initially 
specified  in  the  physical  domain  in  primitive  variables. 

These  are  then  transformed  into  boundary  conditions  in  terms 
of  vorticity  and  stream  function  for  the  computational 
domain. 

Airfoil  Surface.  On  this  surface,  the  no-slip  boundary 
condition  is  enforced,  i.e., 

v  =  0  &  u  =  0  (11) 

It  can  be  seen  that  an  arbitrary  constant  value  of  the 
stream  function  9  satisfies  these  conditions.  In  this  study. 
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the  value  is  chosen  to  be  zero.  Thus,  the  boundary 
conditions  on  the  airfoil  surface  are. 


V  =  0  Vy  =  0  (12) 

As  for  the  vorticity,  because  of  no-slip  at  the  wall, 
vorticity  is  produced  on  this  surface  and  it  is  the  diffusion 
and  subsequent  advection  process,  that  actually  drives  the 
flow  field  [18:140]. 

The  vorticity  at  the  wall  is  computed  using  the 
conditions 

2 An  |-(wn+  V>|  =  -36*  +  46*  -  6*  =  0  (12) 

'  j  ■  i 

V  =  0  (13) 

and  the  recurrence  relations 

6Vj  =  RljfiVj  +  T1j  04) 

6uj  -  R2J6*>J1  +  S2j6wJ.a  +  T2j  05) 

Algebraic  manipulation  of  equations  (12)  to  (15)  results 
in  the  following  relation  for  the  change  in  vorticity  on  the 
airfoil  surface. 


tR13-  4]T12  ♦  S13T22  ♦  T13  -  24W 


[4  -  R 1  ]S1  -  SI  S2 

3  2  3  2 


(16) 


Outer  boundary.  On  this  surface  (Figure  7,  surface  no 
3),  the  free-stream  conditions  due  to  velocity  are  enforced 


u  =  cosa 
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u  =  cosa  and  v  =  sina  (17) 

where  the  non-dimensional  free  stream  velocity  has  unity 
value  and  'a’  is  the  angle  of  attack. 

Applying  the  chain  rule  of  partial  differentiation 

results  in  the  following  condition. 

dtt)  dy  dx  , « _ . 

g*  =  cosa  -g^sina  (18) 

This  is  the  Neumann  form  of  boundary  condition,  and  the 
Dirichlet  form  can  be  obtained  by  integrating  relation  (18), 
which  results  in 

V  =  y  cosa  -xsina+v  (19) 

O 

where  ip  is  the  arbitrary  constant  of  integration  and  in  this 

O 

case  has  already  been  set  to  zero  (airfoil  surface  =  0, 

equation  13). 

Outflow  boundary.  On  this  surface  (Figure  7,  surfaces 
2  and  4),  the  free  stream  conditions  given  by  equations  (17) 
have  also  been  enforced,  which  results  in  the  following 
Neumann  form. 

H  =  §f  eosa  -  If  ,ino  <2°> 

The  corresponding  Dirichlet  condition  is 

ip  =  y  cosa  -  xsina  +  f  (21) 

O 

where  •  =  0 

# 

In  this  study  both  the  Neumann  as  well  as  the  Dirichlet 
type  of  boundary  conditions  have  been  incorporated  in  the 
source  code,  either  of  which  can  be  selected  by  the  user. 
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for  any  surface  (outer  or  outflow).  Based  on  the  author's 
experience  with  the  code,  it  was  been  found  that  the  Neumann 
condition  works  better  on  the  outflow  boundary,  especially 
when  nonsymmetric  flow  is  being  computed.  On  the  outer 
boundary,  the  Dirichlet  condition  works  very  well  and  helps 
accelerate  convergence,  with  greater  stability  than  the 
Neumann  condition.  This  was  found  to  be  especially  true  at 
high  angles  of  attack,  when  application  of  the  Neumann 
condition  can  lead  to  oscillations  and  divergent  behavior. 

It  must  also  be  emphasized  that  each  problem  poses  a 
different  requirement  and  it  is  not  necessary  that  one  type 
of  boundary  condition  will  work  best  under  all  conditions. 
Thus  advantage  of  the  codes  flexibility  can  be  taken  a.  1 
different  boundary  conditions  tested, to  determine  the  optimal 
set  for  the  problem  being  considered. 

Wake  (branch  cut).  Implementation  of  a  boundary 
condition  on  the  branch  cut  poses  a  very  special  problem, 
especially  for  non-symmetr ic  flows,  as  values  of  f  and  u  at 
the  cut  have  to  be  computed.  The  situation  is  further 
complicated  by  the  fact  that  at  the  branch  cut  (C-grid),  a 
discontinuity  of  the  curvilinear  coordinates  exist  (reversal 
of  directions).  Continuity  can  be  maintained,  however,  by 
conceptually  remaining  on  the  same  overlapping  sheet  as  the 
cut  is  crossed.  All  derivatives  thus  do  exist  at  the  cut, 
but  careful  attention  In  difference  formulations  is  necessary 
to  represent  derivatives  correctly  across  the  cut  [25:70]. 
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The  unknowns  y  and  u  are  computed  on  the  branch  cut 
by  solving  the  following  system  of  equations  along 
^-direction  (j=l). 


B>it  C>,.,4  ®>,_,4  E>,4  F>,„=  «;  (21) 

A>.-14  B>,4  C>,.,4  °!6“1.14  E*«u,+  f>141=  R”  (22) 


where  superscript  denotes  a  coefficient  of  the  stream 
function  equation  and  ’o’  denotes  a  coefficient  of  the 
vorticity  equation. 


It  should  be  noted  here  that  the  coefficients  in 
equation  (21)  and  (22)  are  not  computed  in  the  same  manner  as 
for  rest  of  the  domain.  Instead  a  special  finite-difference 
formulation  is  employed  to  enable  handling  of  the  branch  cut 
discontinuity,  as  discussed  in  the  paragraphs  above.  Also, 
the  system  of  equations  is  collated  by  marching  along  the 
constant  rj~line  (j  =  1)  in  the  {-direction.  This  results  in  a 
block  [2x2]  tridiagonal  system,  that  is  solved  by  using  the 
modified  Thomas  algorithm.  Thus,  evaluation  of  *  and  u 
represents  a  two  point  boundary  value  or  elliptic  problem. 
That  must  be  solved  at  each  iteration  level,  based  on  the 
known  boundary  conditions  at  either  end  of  the  branch  cut. 
That  is,  at  the  trailing  edge,  the  boundary  condition  given 
by  equation  (11)  is  enforced  and  at  the  other  end,  the 
boundary  condition  given  by  either  equations  (20)  or  (21)  is 
employed. 
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Initial  conditions 


For  unsteady  flow  solutions  the  proper  initial 
conditions  are  critical.  However,  for  situations  in  which 
only  the  steady  state  solution  is  sought,  the  accuracy  of 
initial  conditions  takes  on  a  slightly  lesser  significance, 
as  the  (theoretically)  steady  state  can  be  achieved  without 
any  regard  to  the  initial  conditions.  But,  due  to  the 
non-linearity  of  the  problem,  it  is  necessary  that  the 
initial  conditions  be  not  too  far  from  the  physics.  Under 
certain  conditions  improper  initial  conditions  may  lead  to 
solutions,  which  may  not  have  any  physical  significance  and 
resemblance  to  the  actual  solution. 

In  this  study  a  second  order  velocity  profile  in  terms 
of  stream  function  is  initialized  in  the  viscous  region  above 
the  airfoil  surface.  In  the  inviscid  region,  the  stream 
function  takes  on  a  distribution  based  upon  the  free  stream 
velocity.  The  vorticity  is  assumed  zero  every  where  in  the 
viscid  as  well  as  the  inviscid  regions. 

Turbulence  modeling 

The  influence  of  turbulence  is  modeled  by  means  of  an 
eddy  viscosity  e,  leading  to  the  total  or  effective  viscosity 

M  +  pc  (23) 

where  jj  is  the  molecular  viscosity. 

The  eddy  viscosity  is  computed  from  the  zero  equation  or 
algebraic  model  of  Baldwin  and  Lomax  [4:2].  The  Baldwin  and 
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Lomax  model  is  patterned  after  the  Cebeci  and  Smith  model 
[101,  which  is  difficult  to  use,  because  of  the  necessity 
of  determining  the  displacement  thickness.  The  Baldwin-Lomax 
model  overcomes  this  short-coming  and  uses  the  distribution 
of  vorticity  to  determine  the  length  scale  in  the  outer 
region  of  the  shear  layer.  Vorticity  is  also  used  in  the 
inner  region,  this  is  beneficial  as  vorticity  is  the  natural 
output  of  the  present  algorithm.  interested  readers  are 
referred  to  [15:6]  for  an  interesting  discussion  on 
shortcomings  of  Baldwin-Lomax  model  and  how  these  can  be 
overcome . 

Details  of  the  Baldwin-Lomax  model  are  discussed  at 
length  in  references  [4:2]  and  [15:6],  therefore  these  will 
not  be  repeated  here. 

The  routine  employed  in  the  computation  of  eddy 
viscosity  is  basically  developed  by  Visbal  [27]  and  has  been 
modified  to  suit  the  present  study.  Also  a  transition  model 
has  been  added  to  account  for  eddy  viscosity  variation,  along 
the  streamwise  direction  in  a  gradual  manner  rather  then  the 
abrupt  change  proposed  by  Baldwin  and  Lomax  [4:2]. 

Transition  model 

Flow  transition  from  laminar  to  turbulent  is  a  very 
complex  phenomenon,  for  which  very  little  understanding 
currently  exists.  The  situation  is  compounded  by  the  complex 
interaction  of  many  factors,  such  as  free  stream  turbulence, 
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surface  roughness,  pressure  and  thermal  gradients,  to  name 
only  a  few.  Hence  prediction  of  transition  relies  very 
heavily  on  experimental  data  and  empirical  formulations.  One 
such  relation  is  due  to  Michel  and  is  known  as  the  Michel’s 
correlation  method  [15:6]. 

Rq  =  1.174^1  +  -L2-.ft£°°  JR""  (24) 

where  subscripts  6  and  x  identify  the  Reynolds  Number  based 
upon  momentum  thickness  and  surface  distance  respectively. 

The  relation  above  is  known  as  the  Michel's  correlation 
formula,  and  is  often  used  to  predict  location  of 
transition.  This  empirical  formula  is  sometimes  coupled  with 
the  location  of  laminar  separation  and  transition  is  assumed 
to  occur  at  either  the  point  of  laminar  separation  or  where 
ever(24)  is  satisfied,  which  ever  comes  first.  It  should  be 
pointed  out  that  Michel's  method  is  not  valid  for  flows  with 
separation  [15:8]. 

In  this  study,  the  transition  location  has  been  picked 
based  on  the  criterion,  presented  below  in  the  order  of 
pref erence . 

(a)  Experimental  data,  if  available. 

(b)  Extrapolation  of  related  experimental  data  reported 

by  Bastedo  [5]  and  Williams  [28]  and  correlation  with 

location  of  laminar  separation. 

(c)  Location  of  laminar  separation. 
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Having  located  the  point  of  transition  the  process  of 
transition  from  laminar  to  turbulent  is  simulated  using  the 
Cebeci-Smith  model  of  transition  [15:6,10],  which  is  given 
by; 

-Q(x  -  ,  j[ 

Jx  * 
tr 

2  3 

Re  U« 

where  Q  =  — T70"0' " ~  1-34  (26) 

where  subscript  tr  designates  values  at  the  beginning  of 
transition  and  i  is  the  transition  factor. 

'tr 

Pressure  and  force  computation 

The  advantage  of  working  with  Vorticity  'V ,  Stream 
Function  ‘V  formulation  is  that  the  number  of  unknown 
variables  reduces,  as  the  pressure  term  is  eliminated  from 
the  primitive  momentum  equation  (see  Appendix  A  for  details). 
Once  <i>  and  ^  have  been  computed,  the  pressure  solution  can 
be  extracted  from  these  known  variables  by  integrating  the 
following  Poisson  equation  for  Pressure  (for  derivations,  see 
Appendix  F). 

Equation  (27)  is  the  Poisson  equation  for  pressure  and 
is  in  general  curvilinear  coordinates.  The  only  unknown 
variable  is  the  pressure  "P"  and  the  right  hand  of  this 
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equation  is  known  from  the  u,  v  solution.  This  equation  is 
descritized  using  central  differences  and  the  following 
equation  is  obtained  at  each  node  point. 


aP+bP+cP  =  d  (28) 

J  J-i  J  J  J  J*i  J 

Solution  along  each  constant  "rj"  line  is  obtained 
implicitly  using  Thomas  algorithm  and  is  advanced  from  one 
level  to  the  next  using  SLOR,  till  convergence  to  a  required 
tolerance  is  achieved. 

Boundary  Conditions  for  Pressure.  The  boundary 
conditions  used  for  the  Poisson  pressure  equation  (27)  are 
mixed,  as  on  the  outer  boundary,  the  Dirichlet  boundary 
condition  P  =  0  is  employed,  the  Neumann  boundary 
condition  0P/6rj  =0  is  enforced  on  the  airfoil  surface  and 
3P/d£  =  0  is  enforced  on  the  outflow  boundary.  On  the 
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branch-cut  (j.i),  the  pressure  is  approximated  by  averaging 
the  value  above  and  below  this  constant  r)  line. 

It  would  be  appropriate  to  point  out  that  at  each 
node  if  P^  C  is  the  solution  obtained  by  integrating 
equation  3.27  (C  is  the  constant  of  integration),  then  the 
level  of  can  be  chosen  by  specifying  the  value  of  P  at 

any  point  in  the  domain.  In  this  study,  the  value  of  P  so 
chosen  is  zero  and  is  specified  at  the  outer  boundary.  This 
approach  facilitates  the  computation  of  the  coefficient  of 
pressure  Cp  ,  which  is  simply  given  by  Cp  =  2 P  .  Thus  the 
variable  P  physically  represents 


P 


P*-  P* 

_ CD 

p!  u*z 


CD 


CD 


where  'V*  represents  dimensional  variables  and  ’V  free 
stream  properties. 


Having  found  the  pressure  and  skin  friction  distribution 
(see  Appendix  Q)  on  the  airfoil  surface,  the  force  and 
moment  coefficients  can  be  computed  using  the  following 
relations,  the  details  on  usage  and  derivation  are  given  in 
Appendix  Q. 


Cl  '  |  As,'(c/  •ir‘(01-a)  -  Cp  cos(e7.-o)j| 

|  As^C^cosCe^-a)  +  Cp  sin(0^-o)j| 

cm  =  "zf  1  A*i(cP  co*ce. ♦  em)  -  c.  ain(er  O.)]) 

c/4  1-1  V.  v  1  JJ 


(29) 
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where  is  the  lift  coefficient,  is  the  dreg  coefficient 

end  C  is  the  moment  coefficient  et  the  querter  chord. 

c/4 

The  other  veriebles  ere  defined  in  Appendix  G 

Convergence  Criterion 

There  ere  three  types  of  converged  solutions  in 
numericel  simuletions,  steady-stete ,  periodic  end  unsteady. 
Even  though  no  solution  is  perfectly  steady,  there  is  always 
some  oscillation.  The  unsteady  behavior  may  be  so  small  as 
to  be  negligible.  When  this  is  the  case  the  solution  can  be 
said  to  have  reached  a  steady-state.  The  steady-state 
solution  is  the  easiest  to  resolve,  and  is  characterized  by 
the  uniform  convergence  of  some  criteria  after  an  initial 
transition  period  [6]. 

The  criterion  for  convergence  for  this  study  were  the 
coefficients  of  lift  (C^)  and  drag  (C^).  When  the  variations 
in  C^  and  C^  are  less  than  0.1*  over  2000  iterations,  then 
the  solution  is  considered  to  have  converged  to  a  steady 
state. 
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IV  Results  and  Discussion 


In  this  chapter  the  Numerical  results  obtained  from  the 
ANS  code  are  presented.  The  verification  of  the  code  and 
grids  is  presented  first,  which  considers  three  cases,  for 
testing  and  evaluating  various  aspects  of  the  code.  This  is 
followed  by  the  results  for  the  Wortman  Airfoil  and 
comparison  with  experimental  data. 

Code  Verification 

Before  the  Wortman  FX  63-137  airfoil  can  be  effectively 
analyzed,  it  is  necessary  to  ensure  that  the  code  and  its 
segments  are  working  correctly  and  are  predicting  the  flow 
field  accurately.  This  has  been  done  by  considering  three 
cases  under  different  flow  conditions  and  comparing  the 
results  so  obtained. 

Case-1;  Laminar  Flow  (Re  =  12,500).  The  first  case 
considered  is  the  symmetric  NACA-12  airfoil  in  laminar  flow 
at  Reynolds  Number  of  12,500. 

The  grid  used  to  study  this  case  was  symmetric  and  had 
outer  boundaries  located  20  chord  lengths  from  the  airfoil 
surface.  The  grid  had  200  symmetrically  placed  grid  points 
on  the  airfoil’s  surface,  spacing  in  (--direction  at  the 
leading  edge  was  0.001  and  at  the  trailing  edge  0.008.  40 

points  were  distributed  on  the  wake  profile,  and  in  the 
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rj-direction  60  points  with  minimum  wall  spacing  of  0.0001. 
Tha  grid  so  obtained  is  shown  in  Figure  11. 

NACA  0012 


-i  a 


ou 


X/C 


10 


20 


Figure  11  Grid  for  NACA-12  airfoil 


The  ANS  code  was  executed  for  this  case  and  convergence 
was  achieved  in  only  9000  iterations,  which  took 
approximately  900  CPU  seconds  to  execute  on  the  Cray  XMP. 
Figure  11  is  a  comparison  of  chord-wise  skin  friction 
distribution  obtained  using  this  code  and  the  results 
obtained  by  Halim  [11],  The  streamline  contours  have  been 
plotted  in  Figure  13a  and  these  can  be  compared  to  contours 
from  reference  [11]  in  Figure  13b,  it  can  be  seen  that 
agreement  so  far  is  very  good.  This  is  further 
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Figure  12  Comparison  of  skin  friction  coefficient 
distribution,  present  with  Ref.  [11]. 


substantiated  by  the  fact  that  laminar  separation  was 
observed  at  x/c  =  0.818,  where  the  location  reported  by 
Halim  [11]  is  x/c  =  0.8178.  Swanson  and  Turkel  [23]  solved 
the  Navier-Stokes  equations  on  a  C-type  grid  for  this  airfoil 
and  Reynolds  number  and  the  obtained  laminar  separation  at 
x/c  =  0.817.  Based  on  the  above  given  results,  it  can  be 
concluded  that  the  code  does  give  excellent  results  for  low 
Reynolds  number  laminar  flow. 

Case-2:  Turbulent  Flow  (Re  =  6.5x10s).  This  test  case 
has  been  designed  to  check  the  accuracy  of  the  turbulence  and 
transition  model,  and  the  ability  of  the  code  to  handle 
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flows  at  high  Reynolds  Number. 

The  geometry  considered  is  the  flat  plate  at  Re=6. 5x10s, 
this  value  of  Re  has  been  chosen  as  turbulent  flow  under 
these  conditions  is  ensured  [20].  The  flat  plate  has  been 
chosen  so  that  the  geometry  of  the  body  does  not  complicate 
the  analysis  of  the  results. 


x/c 

Figure  14  Grid  used  for  Flat  Plate 

To  allow  flat  plate  computation,  the  code  was  slightly 
modified  so  that  a  H-type  grid  may  be  used  (metrics 
singularity  at  the  leading  edge  does  not  allow  use  of  a 
C-type  grid).  The  grid  spacings  used  were  the  same  as  for 
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case-1,  and  Figure  14  is  a  portion  of  the  grid  finally  used 
to  compute  flow. 

To  compare  this  test  case  with  known  and  proven  data, 
the  boundary  layer  solver  (developed  by  Cebeci  &  Smith)  Ref. 
[9]  was  executed  for  Mach  number  of  0.001  and  standard 
atmosphere  free  stream  properties.  This  code  solves  the 
boundary  layer  equations  on  the  infinite  flat  plate, 
advancing  the  solution  parabol ical ly  in  the  etreamwise 
direction,  till  the  required  Reynolds  number  is  reached.  In 
this  case  a  constant  step  size  Ax  =  0.02  was  used  and 
transition  to  turbulent  flow  was  triggered  at  Rx  =  3.2x10^ 

L  20  j . 

The  ANS  code  was  also  executed  based  on  the  above  given 
transition  criterion,  the  differences  being  the  grid 
spacings  and  the  Mach  number,  which  was  zero  as  the  code  is 
limited  to  incompressible  flow  only.  The  results  so  obtained 
are  compared  in  Figures  15  through  17  and  discussed  below. 

Figure  15  compares  the  skin  friction  distribution 
obtained  from  the  two  codes,  and  it  can  be  seen  that  the 
results  do  compare  favorably.  Small  deviation  around  x/c  = 
0.6  can  be  attributed  to  the  differences  in  the  transition 
location  as  it  may  not  locate  exactly  at  the  same  place  due 
to  non-presence  of  an  appropriate  node  point. 

Figure  16  compares  the  velocity  profile  close  to  the 
trailing  edge  (x/c  =  0.98),  with  that  obtained  from  ref. 

[9],  the  two  profiles  agree  very  nicely  and  the  variations, 
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0.015 


Comparison  of  eddy  viscos 
Flat  Plate,  Re  =  6.5x10 


though  very  small  are  within  the  tolerance  (±5%). 

The  third  and  final  plot  for  this  case  is  Figure  17, 
this  graph  compares  the  Eddy  Viscosity  (e)  distribution 
obtained  from  the  two  codes.  It  can  be  seen  that  the 
Original  Baldwin-Lomax  model  [4]  predicts  greater  magnitude 
of  c  in  the  outer  region  and  also  the  distribution  of  the 
Klebanoff  intermi ttency  factor  is  slightly  different,  which 
causes  the  t  profile  to  be  somewhat  different  from  the  one 
obtained  using  Cebeci-Smith  model.  Because  the  Baldwin-Lomax 
model  is  based  on  the  Cebeci-Smith  model  therefore  it  is 
expected  the  later  should  predict  the  same  e  distribution. 

On  order  to  achieve  this  result  the  Baldwin  Lomax  model  was 
modified  with  slightly  different  constants,  given  in  Table 
1.  The  constants  C  is  a  scaling  factor  for  the  outer 

cp 

region  and  is  the  scaling  factor  for  computing  the 

Klebanoff  i ntermi ttency  factor. 

Figure  17  reveals  that  after  modification  the  models  do 
predict  the  same  e  distribution.  it  should  be  noted  that  the 
constants  used  in  this  case  are  applicable  to  only  this 
geometry  and  Reynolds  number.  A  more  general  modification 
however,  is  proposed  by  Visbal  and  Knight  [27]  and  York  [29], 
the  values  of  the  constants  for  the  outer  model  suggested  by 
them  are  given  in  table  1. 

The  constants,  used  in  this  study,  however,  are  the 
ones  given  by  Baldwin-Lomax  ref.  [4],  as  use  of  any  other 
values  would  have  required  a  much  more  in-depth  analysis,  of 
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Table  1  Comparison  of  constants:  modified  and  un-modi fied 
Baldwin-Lomax  Model 


Symbol 

Original 

Present(modif ied) 

Ref. [27] (modified) 

S9 

1.6 

1.33 

1.2 

c 

k  1  •  b 

■ZB 

0.39 

0.646 

the  models  under  different  conditions,  which  is  definitely 
not  within  the  scope  of  this  study. 

Case-3:  Non-Symmetr i c  Turbulent  Flow  (Re  =  200,000). 

The  third  and  final  test  case  is  the  comparison  of  pressure 
distribution  for  Wortman  FX  63-137  airfoil  at  an  angle  of 
attack  of  1°  and  Re  =  200,000. 

The  grid  used  used  to  study  this  case  had  grid  spacing 
and  point  distribution  exactly  the  same  as  for  case-1,  the 
grid  however  so  obtained  was  non-symmetr ic  due  the 
character ist ic  high  camber  of  Wortman  airfoil. 

Figure  18  compares  present  pressure  distribution  with 
data  from  Ref.  [5].  Agreement  between  the  two  is  very  good 
and  variations  are  within  experimental  tolerances. 

Results 

The  results  obtained  for  the  Wortman  FX  63-137  airfoil 
will  now  be  presented,  these  were  obtained  using  the  ANS 
code  developed  in  this  study,  wherein  the  original 
Baldwin-Lomax  model  has  been  used  and  the  surface  of  the 
airfoil  has  been  made  mathematically  smooth,  using  the 
program  "Airfoil"  (Appendix  C). 


Comparison  of  prassura  coafflciant  vs.  x/c 
for  Wortman  airfoil.  Ra  =  ?vin5 


The  grids  used  to  compute  the  different  cases  were  all 
different  from  each  other  primarily  in  the  wake  region.  As 
the  wake  profile  (branch  cut)  had  to  honor  the  trailing  edge 
deflection  and  the  airfoil  angle  of  attack,  it  is  therefore 
not  possible  to  show  here  each  and  every  grid  used  in  this 
study.  However,  Figure  19  is  an  example  of  a  typical  grid 
and  Figure  20  shows  the  grid  details  of  the  leading  as  well 
as  the  trailing  edge.  The  grid  point  distribution  and 
spacings  for  all  grids  were  exactly  the  same  as  defined  for 
case-1  (NACA-12,  Re=12,500). 

Wortman  airfoil,  No  deflections.  Flow  for  this  case  was 
computed  for  6  different  angles  of  attack  ie.  -4,0,4,8,12  and 
16  degrees.  The  constant  stream  function,  velocity  and 
vorticity  contour  plots  for  these  cases  are  attached  in 
Appendix  H  as  Figures  36  through  51.  The  corresponding 
pressure  coefficient  and  skin  friction  distribution  are 
attached  in  Appendix  I  as  Figures  52  through  67.  The 
performance  parameters;  lift  coefficient  (C^),  Drag 

coefficient  (C.)  and  moment  coefficient  (C  )  are  plotted 

a  a 

against  angle  of  attack  in  Figures  21,  22  and  23 
respectively . 

Wortman  airfoil,  5  =  5°,  6  =  0°.  Flow  for  this  case 

■■  1 1  !• 

was  competed  for  5  different  angles  of  attack  ie.  -4, 0,4, 8 
and  12  degrees.  The  constant  stream  function,  velocity  and 
vorticity  contour  plots  for  these  cases  are  attached  in 
Appendix  H  as  Figures  42  through  46.  The  corresponding 
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pressure  coefficient  and  skin  friction  distribution  are 
attached  in  Appendix  I  as  Figures  58  through  61.  The 
performance  parameters;  lift  coefficient  (C^),  Drag 

coefficient  (C,)  and  moment  coefficient  (C  )  are  plotted 

a  « 

against  angle  of  attack  in  Figures  24,  25  and  26 
respectively. 


Wortman  airfoil,  8  =  5°,  8  =  20°.  Flow  for  this  case 

- - 1  m - t« - 

was  computed  for  5  different  angles  of  attack  ie.  -4, 0,4, 8 
and  12  degrees.  The  stream  function,  velocity  and  vorticity 
contour  plots  for  these  cases  are  attached  in  Appendix  H  as 
Figures  46  through  51.  The  corresponding  pressure 
coefficient  and  skin  friction  distribution  are  attached  in 
Appendix  l  as  Figures  62  through  66.  The  performance 
parameters;  lift  coefficient  (C^),  Drag  coefficient  (C^)  and 
moment  coefficient  (C^)  are  plotted  against  angle  of  attack 
in  Figures  27,  28  and  29  respectively. 

Effect  of  Angle  of  Attack 

On  Lift  Coefficient.  The  effect  of  angle  of  attack  on 

the  lift  coefficient  can  be  observed  from  Figures  21,  24  and 

27.  For  all  cases  the  lift  curve  (C^  )  has  a  positive  slope, 

a 

implying  a  lift  coefficient  increase  as  the  angle  of  attack 
is  increased.  It  can  be  also  observed  that  the  lift  curve 
slope  remains  essentially  constant,  however  in  Figure  21  the 
value  of  obtained  for  a=16°  is  less  than  the  value 
predicted  by  the  lift  curve  slope,  this  shows  that  the 
airfoil  is  approaching  the  stall  condition  and  any  further 
increase  in  a  will  trigger  stall  and  would  cause  a 
significant  drop  in  lift.  The  actual  stall  case  could  not  be 
computed  due  to  time  limitation. 

Comparison  with  experimental  results  is  also  made  in 
Figures  21,  24  and  27,  the  present  lift  compare  well  with 
Perry  [16]  for  a  -  -4,  0  and  4  degrees.  For  greater  angles 
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Comparison:  Lift  Coefficiant  vs.  Angle  of  Attack 
Wortman  Airfoil,  6  =0°,  6  =5° 


r.  Drag  Coefficient  vs.  Angle  of  Attack 
rfoll,  6  =0°,  5  =5° 


Moment  Coefficient  vs.  Angle 
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Figure  28  Comparison:  Drag  Coefficient  vs.  Angle  of  Attack 
Wortman  Airfoil,  5  =20°.  6  =5° 


□  Ref  ( 

©PRES 
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Worti 


of  attack  the  disparity  is  larger  and  is  more  than  the 
acceptable  limits.  Figure  21  also  compares  the  present  data 
with  data  reported  by  Williams  [28]  for  the  no  deflections 
case  and  various  Reynolds  numbers.  It  can  be  seen  that  this 
data  is  also  significantly  different  from  Ref  [16].  For  Ref. 
[b28]  the  lift  curve  slope  is  greater,  however  in  the 
present  case,  the  lift  curve  slope  is  lower  than  Perry  [16]. 
Another  point  to  note  is  that  airfoil  stalls  at  quite 
different  values  of  a  for  Ref.  [16]  and  Ref  [28]  (The  present 
results  reflect  the  uncertain  data  scattering  of  experimental 
measurements,  at  low  Reynolds  conditions). 

On  Drag  Coefficient.  The  drag  coefficient  has  been 

plotted  as  a  function  of  angle  of  attack  in  Figures  22,  25 

and  28,  The  drag  distribution  is  not  linear  ,  but  is  closer 

to  parabolic.  Because  of  the  limited  o  range  considered  the 

variation  of  C.  beyond  the  tested  values  cannot  be  predicted. 

a. 

In  comparison  to  the  other  cases,  the  data  set  corresponding 
to  leading  edge  deflection  of  5°  is  the  closest  to  the  data 
reported  by  Perry.  For  all  cases  the  variance  between 
present  and  Perry  is  greatest  at  higher  angles  of  attack,  in 
comparison  to  values  of  o  close  to  zero,  where  the  agreement 
is  very  good.  Figure  22  also  compares  variation  with  Ref. 
[28]  the  values  reported  are  an  order  of  magnitude  lower  than 
those  reported  by  Perrv.  In  fact  the  magnitude  difference  in 
measurements  between  [16]  and  [28]  is  similar  to  the 
disparity  between  present  calculations  and  data  of  Ref. [16]. 
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On  Moment  Coefficient.  Positive  pitching  moment  is 


defined  in  this  study  as  one,  which  causes  a  nose  down 
motion  of  the  airfoil.  In  all  the  cases  considered  the 
variation  of  moment  coefficient  with  a  is  almost  linear. 

This  can  be  clearly  seen  in  Figures  23,  26  and  29,  along 
with  the  fact  that  agreement  with  experimental  data  is  pretty 
good  within  the  a  =  ±  4  degree  range.  For  higher  values  of 
angle  of  attack,  that  is  8  and  12  degrees,  no  agreement 
between  data  and  calculations  is  reached.  Similar  to  the  the 
drag  and  lift  coefficients,  the  case  with  leading  edge 
deflected  of  5  degrees  is  closest  to  experimental  data,  in 
comparison  with  the  other  two  cases.  The  attention  should 
also  be  focused  on  the  fact  that  the  data  from  Ref.  [28]  is 
quite  different  than  Ref.  [16].  The  significant  discrepancy 
in  experimental  measurements  and  between  numerical  simulation 
suggest  a  continuous  research  effort  in  this  area  are  still 
required. 

On  Pressure  Peak.  Close  examination  of  the  Cp  plots 
Figures  52  through  67  shows  that  at  smaller  angles  of  attack, 
it  is  very  hard  to  pick  out  the  location  at  which  minimum 
pressure  occurs  as  the  Cp  variation  is  quite  flat.  However, 
as  the  angle  of  attack  is  increased  the  point  of  min. 
pressure  becomes  more  and  more  distinct  and  moves  upstream, 
towards  the  leading  edge. 
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Table  2  Comparison  of  separation  and  reattachment  points 
Wortmen  FX  63-137  airfoil,  6  =5  =0° 

te  le 


Present 

0.75 

-» 

X 

0.60 

-» 

X 

0.65 

-» 

0.72 

0.95 

X 

0.49 

-> 

X 

0.43 

-♦ 

0.51 

0.60 

-» 

X 

0.17 

-» 

0.36 

0.38 

-♦ 

0.47 

0.48 

X 

(a)  Lower  surface 


Bastedo 
Re=2x 10s 
Ref.  [5] 


0.59 


0.48 


0.31 


0.82 


0.72 


0.60 


Willi ams 
Re=7x10s 
Ref.  [28] 

0.63 

-> 

0.78 

0.52 

-* 

0.70 

0.47 

-* 

0.64 

0.40 

-* 

X 

0.11 

-¥ 

0.16 

0.75 

-» 

X 

0.05 

-¥ 

0.12 

0.50 

-* 

X 

(b)  Lower  surface 


Present 

Bastedo 
Re=2x 10s 
Ref.  [5] 

Willi ams 
Re=7x105 
Ref.  [28] 

00  *- 
(O  <J> 

•  • 

o  o 

t  T 

CO  CM 
00 

o  o 

0.05  -*  0.08 

0.05  -♦  0.08 

0.32  -»  0.81 

0.38  -♦  0.52 

0.40  -»  0.49 

0.66  -»  0.67 
0.95  -*  0.99 

00 

CO 

• 

o 

t 

CO 

• 

o 

0.48  -»  0.52 

0.53  -»  0.76 

0.60  -»  0.78 

0.62  ->  0.68 

NIL 

NIL 

NIL 

NIL 

NIL 

NIL 
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The  same  deductions  can  be  arrived  at  by  observing  the 
constant  velocity  contour  plots,  given  in  Figures  36  through 
51,  where  it  can  be  seen  that  as  the  angle  of  attack  is 
increased,  the  point  of  minimum  pressure,  which  is  also  the 
point  of  maximum  velocity  moves  towards  the  leading  edge. 

A  point  to  note  in  Figures  59  through  67  is  that  the 
pressure  coefficient  shows  a  small  spike  of  same  magnitude 
and  is  located  at  the  same  place  (x/c  =  0.13)  for  all  cases, 
with  slat  deflection.  This  spike  has  been  generated  due  the 
slat  deflection,  and  can  be  removed  by  smoothing  the  metric 
discontinuity  at  this  point  on  the  airfoil. 

On  Laminar  Separation.  The  point  of  laminar  flow 
separation  is  defined  by  the  location  of  the  vanishing  skin 
friction  coefficient.  The  relation  of  the  point  of 
separation  with  angle  of  attack  can  be  observed  from  the 
stream  line  plots  given  in  Figures  36  through  51.  The  exact 
location  of  the  separation  point  can  be  found  from  the  skin 
friction  distribution,  which  is  plotted  as  a  function  of 
chord  in  Figures  52  through  67.  Note  that  at  the  point  of 
separation  the  values  of  the  stream  function  as  well  as  the 
vorticity  are  zero. 

From  the  stream  function  contours  it  can  be  very  easily 
observed  that  as  angle  of  attack  is  increased,  the  point  of 
separation  on  the  upper  surface  moves  upstream  towards  the 
leading  edge  and  the  separation  bubble  becomes  bigger  in 
length  and  height.  On  the  lower  surface  however,  the  point 


of  separation  moves  downstream  and  the  separation  bubble 
becomes  smaller  and  smaller  in  length  till  it  completely 
vanishes.  The  fore  mentioned  phenomenon  can  be  explained  by 
the  fact  that  on  the  upper  surface  the  adverse  pressure 
gradient  becomes  stronger  as  a  is  increased,  which  triggers 
earlier  separation.  On  the  lower  surface,  the  opposite 
takes  place  as  flow  on  this  side  has  greater  momentum  which 
helps  to  overcome  the  adverse  gradient  causing  the  point  of 
separation  to  move  downstream.  This  trend  is  most  evident 
from  the  data  presented  in  Tables  2  through  4  .  These  tables 
compare  the  values  of  separation  and  reattachment  for  the 
present  study,  with  flow  visualization  results  of  previous 
experiments,  wherever  available  (a  x  identifies  no 
reattachment) . 


Table  3  Comparison  of  separation  and  reattachment  points 
Wortman  FX  63-137  airfoil  6  =5°,  6  =0° 
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0 

0.40  -»  0.59 

0.58  -»  0.99 
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Table  4  Comparison  of  separation  and  reattachment  points 


Wortman  FX  63-137  airfoil  6  =5°,  6  =20° 

1  •  !• 


o 

a 

Present 

Perry 

Re=1x10 

Ref.  [16] 

Lower  surf. 

Upper  surf. 

Upper  surf. 

-4 

0.01  ->  0.10 

0.53  -»  0.91 

0 

0.37  -+  0.75 
0.21  -»  0.73 

0.57  ->  0.90 

4 

0.36  -4  0.74 

0.49  -»  0.84 

0.60  -»  0.70 

8 

0.36  -4  0.82 

0.27  -+  0.51 
0 . 66  0*96 

0.35  —4  0.60 
0.75  -4  x 

12 

0.44  ->  0.77 

0.41  -4  0.78 
0.83  -4  0.91 

0.20  -4  0.30 
0.75  -4  x 

Effect  of  Slat  Deflection 

The  slat  constitutes  13%  of  the  chord  and  in  this  case 
has  been  deflected  by  an  angle  of  5  degrees.  Downward 
deflection  of  slat  is  defined  positive. 

On  lift  coefficient.  The  effect  can  be  best  evaluated 
by  observing  Figure  30,  which  is  the  comparison  of  lift 
curve  slopes  for  the  three  different  cases.  Comparing  the 
case  of  no  slat  deflection  with  the  case  of  5  degree 
deflection,  it  can  be  seen  that  at  small  angles  of  attack 
lift  coefficient  for  the  Wortman  Airfoil  with  slat  deflection 
is  slightly  less  then  compared  to  the  no  deflection  case. 
However,  at  higher  angles  of  attack,  the  opposite  is  true. 
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The  above  observations  have  also  been  reported  by  Perry 
[16],  who  points  out  that  at  deflection  of  slats  increases 
the  value  of  the  critical  angle  (angle  at  which  an  abrupt 
change  of  parameters  takes  place;  stall). 

These  observations  can  be  best  explained  by  the 
following;  at  lower  angles  of  attack  and  slat  deflection  , 
the  effective  angle  of  attack  seen  by  the  airfoil  leading 
edge  decreases,  causing  reduction  of  the  pressure  peak  and 
it's  movement  downstream.  The  net  results  is  a  reduction  in 
lift.  At  higher  angles  of  attack,  however,  downstream 
movement  of  the  point  of  minimum  pressure  delays  laminar 
separation  and  hence  stall,  this  improves  airfoil 
performance  and  increases  the  critical  angle. 


Figure  30  Effect  of  slat  and  flap  deflections  on  lift 
coefficient 
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On  Drag  Coefficient.  The  effect  of  slat  deflection  on 
the  drag  coefficient  is  made  obvious  by  comparing  the  case  of 
no  deflection,  with  the  case  of  5  degree  slat  deflection  in 
Figure  31.  Almost  no  change  in  c^  occurs  due  to  the  slat 
deflection.  Perry  [16]  also  reports  a  similar  finding,  for 
the  angle  of  attack  range  considered.  The  effect  of  slat 
deflections  of  10  degrees  and  greater  is  felt  only  when  the 
angle  of  attack  is  greater  than  16  degrees  and  less  than  -8 
degrees . 


Figure  31  Effect  of  slat  and  flap  deflections  on  drag 
coefficient 
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On  Moment  Coefficient.  For  the  leading  edge  deflection 
of  5  degrees,  the  effect  on  moment  coefficient  is 
negligible.  This  can  be  observed  from  Figure  32,  where  the 
moment  coefficient  for  the  deflected  case  is  compared  to  the 
un-deflected  case.  Perry  [16]  also  reports  a  similar  finding 
for  the  deflection  angle  and  angle  of  attack  range,  being 
considered  here. 


Figure  32  Effect  of  slat  and  flap  deflections  on  moment 
coefficient 

On  Pressure  Peak.  Comparison  of  the  pressure 
coefficient  plots  for  the  undeflected  case  in  Figures  52 
through  54  with  the  slat  deflected  case  in  Figures  55  through 


75 


57  shows  very  small  localized  changes  in  the  leading  edge 
zone.  The  observation  is  anticipated,  that  is,  slat 
deflection  causes  the  point  of  minimum  pressure  to  move  down 
stream  and  reduces  the  pressure  peak.  For  all  the  cases  with 
slat  deflection,  a  spike  in  Cp  distribution  to  appears  at  the 
slat  rotation  point.  This  anomaly  is  induced  by  the  metric 
discontinuity  at  that  point  on  the  upper  surface.  This  minor 
discrepancy  can  be  eliminated  by  smoothing  the  upper  surface 
at  that  point. 

On  Separation  Point.  Effect  of  the  slat  deflection  on 
the  separation  points  can  be  evaluated  by  comparing  the 
separation  points  given  in  Table  3  for  slat  deflection  with 
separation  point  data  given  in  Table  2  for  the  un-deflected 
case.  There  is  negligible  change  in  the  location  of  the 
separation  point.  This,  however,  would  not  be  true  for  large 
slat  deflections. 


Effect  of  Flap  Deflection 

The  Flap  constitutes  25*  of  the  chord  and  in  this  case 
has  been  deflected  by  an  angle  of  20  degrees.  The  downward 
deflection  of  the  flap  has  been  defined  positive. 

On  Lift  Coefficient.  There  is  a  significant  increase  in 
the  lift  coefficient  due  to  flap  deflection,  depicted  in 
Figure  30.  The  complete  lift  curve  slope  has  been  shifted 
upwards,  this  is  the  usual  character istic  of  flap 
deflection.  Due  to  non-linearity  and  viscous  effects,  there 
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has  been  a  very  slight  change  in  the  lift  curve  slope  also. 
The  flap  deflection  has  increased  the  slope  by  a  very  small 
value.  It  can  be  predicted  that  further  increase  in  flap 
deflection  angle  will  cause  more  shift  of  the  curve  in  the 
same  direction. 

This  increase  in  the  lift  coefficient  can  be  explained 
by  the  constant  stream  function,  vorticity  contour  plots 
(Figure  47  -  51)  and  the  coefficient  of  pressure  graphs 
(Figure  63  -  76).  Due  to  flap  deflection,  the  pressure 
distribution  has  been  drastically  changed.  In  the  flap 
region  a  large  suction  exists  on  the  upper  surface,  on  the 
lower  surface  positive  pressure  has  been  also  increased. 

This  combined  effect  has  been  created  due  to  the  counter 
rotating  flow  occupying  the  void  left  by  the  deflected  flap. 
This  counter  rotating  flow  can  be  observed  in  Figure  33, 
which  shows  the  velocity  vector  plots  for  two  of  the  flap 
deflected  cases.  This  counter  rotating  vorticity  strengthens 
the  circulation  around  the  airfoil,  which  increases  the 
airfoi 1  lift. 

On  Drag  Coefficient.  The  change  in  drag  coefficient  due 
to  flap  deflection  is  presented  in  Figure  31.  The  complete 
curve  has  been  shifted  up  by  a  certain  constant  value.  The 
viscous  effects  and  the  non-linearities  associated  with  this 
flow  field  have  caused  a  slight  change  in  the  drag 
distribution . 

Most  of  the  increase  in  drag  is  due  to  the  pressure 
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(b)  a  =  12° 

Figure  33  Velocity  vector  plots,  Wortman  airfoil 
6  =5°,  6  =20° 

l  •  • 
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differential  (explained  above),  acting  on  the  flap  surface, 
with  flap  deflected,  the  wetted  area  (effective 
crossect ional  area  seen  by  the  free  stream)  has  also 
increased,  which  also  contributes  to  an  increase  in  drag. 

On  Moment  Coefficient.  The  same  phenomenon  occurs  for 
the  moment  coefficient  as  has  been  observed  for  the  above  two 
cases.  Figure  32  shows  that  the  moment  coefficient  curve  has 
been  shifted  by  a  certain  factor,  with  a  small  change  in 
slope. 

The  increase  of  moment  coefficient  can  be  explained 
again  as  being  the  result  of  the  large  pressure  differential 
acting  on  the  airfoil  surface,  specifically  in  the  flap 
region  (the  complete  pressure  distribution  around  the  airfoil 
changes,  to  accommodate  the  increase  in  circulation).  Due  to 
the  long  moment  arm,  even  a  small  change  in  the  pressure 
distribution  in  the  trailing  edge  zone  result  in  significant 
changes  in  the  moment  coefficient. 

On  pressure  peak.  Comparing  the  pressure  distributions 
of  the  flap  deflected  (Figures  64  -  67),  with  the 
un-deflected  case  (Figures  58  -  63)  for  all  angles  of  attack, 
the  trend  is  moving  the  point  of  minimum  pressure  towards  the 
leading  edge.  A  considerable  increase  in  adverse  pressure 
gradient  on  the  upper  surface  and  corresponding  pressure  rise 
on  the  lower  surface  are  also  detected. 

On  Separation  Point.  Due  to  the  forward  movement  of  the 
point  of  minimum  pressure  and  the  increase  in  adverse 


pressure  gradient  acting  on  the  upper  surface,  the  point  of 
separation  occurs  much  closer  to  the  leading  edge.  This  can 
be  observed  from  the  constant  stream  function  contours 
(Figures  47  -  51)  and  comparison  of  separation  data  given  in 
Table  4  with  separation  data  given  in  Tables  2  and  3.  From 
the  stream  function  and  vorticity  contours,  note  that  as  the 
angle  of  attack  increases  with  flap  deflection,  the  region 
of  separated  flow  becomes  larger  with  point  of  separation 
further  upstream.  As  flap  deflection  or  angle  of  attack  is 
increased,  the  separation  will  become  larger  and  will 
subsequently  cause  stall  of  the  airfoil. 
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Conclusions  And  Recommendations 


The  results  obtained  using  the  ANS  code  developed  in 
this  study  compared  favorably  with  experimental  data  for 
small  angles  of  attack.  For  larger  angles,  the  agreement 
needs  improvement.  However  based  on  the  codes  performance, 
it  can  be  said  that  the  solutions  obtained  are  still  good 
approximation  of  the  physical  problem.  It  has  also  been 
observed  that  experimental  data  also  varies  from  source  to 
source  and  the  discrepancies  are  comparable  in  magnitude,  in 
comparison  with  numerical  results  obtained  using  the  ANS 
equations . 

The  results  indicate  that  deflection  of  the  flaps  causes 
earlier  separation  of  the  boundary  layer  and  slat  deflection 
tends  to  negate  this  effect.  Because  only  a  limited  number 
of  cases  were  computed,  no  optimization  studies  could  be 
carried  out. 

Convergence  character istics  of  the  solution  is 
reasonable  for  the  Reynolds  number  considered  here,  however, 
the  code  does  have  a  Reynolds  number  limitation.  No  studies 
were  carried  out  to  find  that  limiting  value,  as  it  would 
vary  from  problem  to  problem.  The  code  also  exhibited  some 
sensitivity  to  the  initial  conditions,  improper 
initialization  of  the  flow  field  causes  delayed  convergence 
in  limit  cases.  This  is  because  of  the  fact  that  the  stream 
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function  and  vorticity  equations  are  being  solved  in  a 
coupled  manner.  The  problem  can  be  eliminated  if  the  stream 
function  equation  is  solved  using  direct  methods  at  each  time 
or  iteration  level.  The  recommended  procedure  may  alleviate 
the  sensitivity  to  initial  conditions  for  convergence.  It  is 
also  recommended  that  the  ADI  scheme  may  be  used  in 
conjunction  with  artificial  damping,  which  will  ultimately 
allow  reduction  in  the  Reynolds  Number  limitation  and  will 
also  improve  rate  of  convergence.  It  is  also  recommended 
that  Baldwin-Lomax  model  be  evaluated  together  with  the 
Cebeci  Smith  model  for  turbulence  for  the  investigated  cases. 
It  is  felt  that  the  large  disparity  at  higher  angles  of 
attack  is  induced  by  the  turbulence  model  to  effectively 
handle  the  massive  separations  encountered  in  this  study. 

Finally  it  can  be  concluded  that  the  ANS  code  holds  good 
promise  and  refinements  in  the  designated  areas  will  allow 
its  development  into  a  good  design  tool. 
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APPENDIX  A  :  Derivation  of  stream  function  and  vorticity 


transport  equations 

The  non-dimens ional i zed  Navier  Stokes  equations  are 


du  .  3v 
3x  5y 


f¥  *  ufSi  +  vfH 

at  3x  dy 


dv  ,  dv  .  dv 

WE  +  ^  +  va* 


ap  .  _L  af  ai/j  +  a  fain 
a*  +  aS^J  +  ay^J 

dp  .  i  a  f  av\  .  a  f  evi 
5?  +  “r5  ax^^SxJ  +  a?^ 


(30) 

(31) 

(32) 


Differentiating  equation  (31)  with  respect  to  y  and 
subtracting  equation  (32)  differentiated  with  respect  to  x 
from  the  first  elliminates  the  pressure  term  and  the 
resulting  equation  reduces  to  the  following. 


|^(Eqn.  32)  -  ^(Eqn.  31)  : 


a  (dv  aul  „a  fav  au)  a  fav  aui 
WE{Wi  Wy)  +  “^[53?  "  Wy)  vWy[Wx  "  5yJ 

.  i  faz_f..av\  a2  f  au\  A  a2  f  av|  a2  f  au\l 
■  “Re  |^x2[M5yJ  '  3y3x['i3xJ  +  5y5x[M5yJ  "  Wy2[^y)j 


(33) 


applying  the  defination  of  vorticity  and  stream  function  to 
the  left  hand  side  (LHS)  of  above  equation  results  in  the 
following  form; 
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* 


(34) 


,  _  do  .  dw  dfa>  dtt  du 

LHS  "  +  5y  33?  '  33?  3y 

whsre  u  =  §5  -  gp  (35) 

«nd  u  =  §t  »  =  -  I*  (36) 

the  right  hand  side  can  be  expanded  term  by  term  as  follows 
lx*  ('■§;)  5  “  yxxx  +  "x  vxx  +  »,x  vx  +  *V  vxx 

■  ^  —  I  -  LI  U  +JJU  +  U  U+UU 

aydx(^3xj  H  xxy  Hx  xy  *xy  x  My  xx 
=  "  vyyx  *  “x  vyy  +  •'xy  vy  *  "y  vxy 

I A"*})  =  “  Vy  +  "y  uyy  +  V  uy  *  "y  uyy 


substituting  the  above  expansions  into  RHS  of  equation  (33) 
results  in  the  following  simplified  form 

RHS  I-S«(RJxx  +  *"yy  +  "x“x  +  »‘y“y  +  "x  'xx*  "xx  vx'  Ryuxx 

+  "x  vyy  tiyuyy  ^ yy  uy  Rxyux  +  ^xy  vy} 

making  use  of  the  definations  of  continuity  (30)  and 
vorticity  (35)  the  above  expression  can  be  recast  as  follows. 


RKS  =  ♦  2  »x»x  * 

U  fa)  +  Ufa) 

Hxx  H yy 

+  2 

Ufa)  +  U  ( 

Hy  y  Hyy 

+ 

2  u 

v  + 

u 

u  -  u 

w  y 

Hxy 

yy 

Hxx 

yy  x 

it  can  be  seen  that  the  above  expression  is  the  expanded  form 
of  the  following  relation 
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[fcH]  *  ly[§?H]  "  2 5x5/ 


2  _  2 
JL  d  £ 
6y  3x3y 


+  a*n  a2»  +  a2/i  a2y> 
3x9x  ay3y  dydy  dxdx 


Thus  the  vorticity  transport  equation  is: 

3w  .  dtp  6u  dtp  du  _  i  /  a  pa  r  .  a  ("a  f  H 
aS  *  a?  5*  3?ay--S5|  aslMHj  a/b/MJ 

,  a 2m  a2*  .  a2i*  a2v  .  a2n  a2* 

dxdy  dxdy  dxdx  dydy  dydy  dxdx 


(37) 


Substituting  the  defination  of  stream  function  into  the 
continuity  equation  causes  it  to  vanish,  therefore  the  stream 
finction  equation  is  obtained  by  subtituting  (36)  into  (35), 
which  results  in  the  following  form. 


+  JOl  =  „ 

dxdx  dydy 


(38) 


Summar i z ing 


Vorticity  transport  equation 


aw  .  dw  a«  dtp  aw  _  i  /  a  fa  f  11  .  a  fa  f  1 

d~E  +  dy<5x”dxFy”  ~ Re  |  +  SylMH 


_  o  dZii  dZtp  .  d~u  d~tp  x  d~n  d~tp  , 

23S5y  3*3y  +  3S3x5]73y  +  373y  (39) 


a2u  a2 


a2u  a2* 


Stream  function  equation 


a  y>  .  6u>  _ 

dxdx  dydy  ~  u 
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APPENDIX  B 


Co-ordinate  Transformation  Relations 
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where  f  is  any  arbitrary  function 
and 

J  =  XE  Yr,  '  Xr,  yE 

2  .  2 

°  "  yn  xn 

0  =  y  y_  +  x  x_ 
H  5  rj  5 

2  .  2 

»  =  yE  +  XE 


2  2 
•  =  yn  ’ 


P  =  y»  yE  '  X1  XE 
'  2  2 
»  =  ye  -  xs 


where  subscripts  denote  dif ferent iation  w.r.t.  variable 
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APENDIX  C  :  Fortan  listing  of  program  "AIRFOIL" 


PROGRAM  AIRFOIL 

- PROGRAM  TO  DEFLECT  THE  COORDINATES  OF  THE 

AIRFOIL’S  LEADING  AND  TRAILING  EDGE,  AND  SMOOTH  THE  DATA 

DIMENSION  X(301),Y(301) , XU (200) , YU (200) , XL (200) , YL(200) 
CHARACTER* 15  INPUT 
C 

WR ITE(* ,*) '  NAME  OF  FILE  CONTAINING  AIRFOIL  DATA' 


READ(*, 1) INPUT 
1  FORMAT (A) 

OPEN ( 1 , F I LE= I NPUT ) 

C - READ  AIRFOIL  CONFIGURATION  DATA 


WR I TE (*,*)'  SLAT  LENGTH  (CHORD  RATIO)' 
READ(* , * )SLAT 

WRITE(*,*)'  FLAP  LENGTH  (CHORD  RATIO)’ 
READ(* ,*)FLAP 

WR I TE (* , * ) '  SLAT  DEFLECTION  (DEGREES)' 
READ(* ,*)SDA 

WR I TE(* ,*) '  FLAP  DEFLECTION  (DEGREES)' 
READ(* ,*)FDA 
Pl=3. 141592654 
SDA=SDA*P 1/180. 

FDA=FDA*P 1/180. 

AF  =  1.0-FLAP 

C - READ  AIRFOIL  SURFACE 

KOUNT =0 

5  CONTINUE 

READ( 1 , * , END= 10 )X( KOUNT  + 1 ) ,Y(KOUNT+1 ) 

KOUNT=KOUNT+1 

GOTO  5 

10  CONTINUE 

CLOSE ( 1 ) 
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NBODY=KOUNT 

WR I TE(* , *) '  NO  OF  POINTS  READ  IN,  NBODY= ' , NBODY 

c - SET  LEADING  EDGE  POINTER 

XMIN=1000 

DO  20  I  =  1, NBODY 

I F(X( I )  .LT.  XM I N )THEN 
XMIN  =  X( I ) 

JLE  =  I 
END  IF 


20  CONT I NUE 

C - LOWER  SURFACE  POINTER 

KL=0 

DO  30  I  =  JLE,  1,-1 


XL(KL+1)=X( I ) 
YL(KL+1 )=Y( I ) 
KL=KL+ 1 


30  CONT I NUE 

C - UPPER  SURFACE  POINTER 

KU=0 

DO  40  I  =  JLE, NBODY 


XU(KU+ 1 )=X( I ) 
YU(KU+1)=Y( I ) 
KU=KU+ 1 


40  CONTINUE 
C 

C================================SLAT================= 

C 

C - SET  SLAT  POINTERS - 

C - (UPPER  SURFACE) 

I F ( (SLAT  .EQ.  0.0)  .OR.  (SDA  .EQ.  0.0))GOTO  5000 
DO  50  I  =  1 , KU 

1 F (  XU(I)  .EQ.  SLAT)  TSU  =  1.0 
1 F (  XU( I )  .LT.  SLAT)  THEN 
JSU=I 
END  IF 
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50  CONTINUE 

c - (LOWER  SURFACE) 

DO  60  I  =  1 ,KL 


I F (  XL( I )  .EQ.  SLAT)  TSL  =  1.0 
I F (  XL ( I  )  .LT.  SLAT)  THEN 
JSL=  I 
END  IF 

60  CONTINUE 

c - GENERATE  POINTS  CORRESPONDING  TO  SLAT 

YUSLAT=YU( JSU)+(YU( JSU+1 ) -YU (JSU) )*(SLAT-XU( JSU) ) 
1  /(XU(JSU+1)-XU(JSU)) 

YLSLAT=YL ( JSL ) + ( YL ( JSL+ l)-YL(JSL) )* ( SLAT-XL ( JSL ) ) 
1  / ( XL ( JSL+ 1 ) -X_ ( JSL  )  ) 

c - (UPPER  SURFACE) 

c - CREATE  CAVITY  IN  ARRAYS 

I F (  TSU  .EQ.  1.0)  THEN 
JSU=JSU+1 
ELSE 


DO  70  I  =  KU , JSU+ 1,-1 
YU( l+1)=YU( I ) 

XU( l+1)=XU( I ) 

70  CONTINUE 

C - INSERT  INTO  ARRAY 

YU( JSU+1 )=YUSLAT 
XU( JSU+1 )=SLAT 

C - RESET  POINTERS 

KU=KU+1 
JSU= JSU+1 


END  IF 


c - (LOWER  SURFACE) 

c - CREATE  CAVITY  IN  ARRAYS 


I F (  TSL  .EQ.  1.0)  THEN 
JSL= JSL+ 1 
ELSE 

DO  80  I  =  KL, JSL+1 ,-1 
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YL( l+1)=YL( I ) 

XL ( 1+1 )=XL ( I ) 

80  CONTINUE 

C - INSERT  INTO  ARRAY 

YL ( JSL+ 1 ) =YLSLAT 
XL ( JSL+ 1 ) =SLAT 

C - RESET  POINTERS 

KL=KL+ 1 
JSL= JSL+ 1 


END  IF 


C 

c - ROTATE  SLAT  ACCORDING  TO  DEFLECTION  ANGLE 

C  APPLY  COORDINATE  TRANSFORMATION 

c - UPPER  SURFACE 


XS=SLAT 

YS=YLSLAT  + ( YUSLAT-YLSLAT )/2 . 0 
SHYP=SQRT(XS*XS+YS*YS ) 

SANG=ATAN2(YS,XS) 

DO  130  I  =  1,JSU 

CALL  SROTAT ( XS , YS , SHYP , SANG , SDA , XU( I )  ,YU( I ) ,XX,YY) 
XU( I )=XX 
YU ( I )=YY 


130  CONTINUE 

c - LOWER  SURFACE 

DO  140  I  =  1 , JSL 


CALL  SROTAT(XS,YS, SHYP, SANG, SDA, XL ( I ) , YL ( I ) , XX , YY ) 
XL ( I )=XX 
YL ( I ) = YY 


140  CONTINUE 

I F ( SDA  .GT.  0.0)  THEN 

C - POSITIVE  SLAT  DEFLECTION 

C - - - REMOVE  OVERLAPPING  POINTS  FROM  LOWER 

C  SURFACE 

C - FIND  HOW  MANY  POINTS  OVERLAP 

I  FO=0 
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DO  110  I  =  1 , JSL 

I F (XL ( I )  .GT.  SLAT)  THEN 
IFO= IFO+1 
END  IF 

110  CONTINUE 

c - REMOVE  IFO  NUMBER  OF  POINTS  FROM  LOWER 

C  SURFACE 

JTEMP= JSL- IFO+1 
DO  120  I  =  JTEMP , KL- IFO 
XL ( I )=XL( I+IFO) 

YL( I )=YL( I+IFO) 

120  CONTINUE 

C - RESET  POINTERS 

KL=KL- IFO 
JSL=JSL- IFO 

c - REMOVE  POINTS  FROM  UPPER  SURFACE  FOR 

C  CURVATURE  CONTINUITY 


IF  ( ABS(SDA)  .GT.  5*PI/180.0)  THEN 
DO  181  I  =  JSU- 1 , KU- 1 
XU( I )=XU( 1  +  1  ) 

YU( I )=YU( 1  +  1  ) 

181  CONTINUE 


KU=KU-1 
END  IF 
ELSE 

c - NEGATIVE  SLAT  DEFLECTION 

c - REMOVE  OVERLAPPING  POINTS  FROM  UPPER 

C  SURFACE 

C - FIND  HOW  MANY  POINTS  OVERLAP 

I  F0=0 

DO  111  I  =  1,  JSU 


I F (XU( I )  .GT.  SLAT)  THEN 
!FO= IFO+1 
END  IF 

1 1 1  CONT I NUE 
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■REMOVE  IFO  NUMBER  OF  POINTS  FROM  UPPER 
SURFACE 


JTEMP=JSL- IFO+1 
DO  121  I  =  JTEMP , KU- I FO 
XU( I )=XU( I+IFO) 

YU( I )=YU( I+IFO) 


CONT I NUE 


■RESET  POINTERS 


KU=KU- I FO 
JSU= JSU- I FO 


- REMOVE  POINTS  FROM  LOWER  SURFACE  FOR 

CURVATURE  CONTINUITY 
IF  (ABS(SDA)  .GT.  5*PI/180.0)  THEN 
DO  182  I  =  JSL-1,KL-1 
XL ( I )=XL( 1+1) 

YL( I )=YL( 1+1) 

CONT I NUE 
KL=KL- 1 
END  IF 
END  IF 


5000  CONTINUE 


C===== ===========-======== =========FLAP: 


■SET  FLAP  POINTERS- 


- (LOWER  SURFACE) 

IF((  FLAP  .EQ.  0.0)  .OR.  (FDA  .EQ.  0.0))GOTO  6000 
DO  51  I  =  1  ,  KL 

I F (  XL( I )  .EQ.  AF)  TFL  =  1.0 
I F (  XL ( I )  .LT.  AF)  THEN 
JFL=  I 
END  IF 
CONT I NUE 


- (UPPER  SURFACE) 

DO  61  I  =  1 , KU 


i 


I F (  XU( I)  .EQ.  AF)  TFU  =  1.0 
I F (  XU( I )  .LT.  AF)  THEN 
JFU=  I 
END  IF 

61  CONTINUE 

C - GENERATE  POINT  CORRESPONDING  TO  FLAPS 

YUFLAP=YU(JFU)+(YU(JFU+1)-YU(JFU))*(AF-XU(JFU)) 

1  / (XU( JFU+1 )-XU( JFU) ) 

YLFLAP=YL(JFL)+(YL(JFL+1)-YL(JFL))*(AF-XL(JFL)) 

1  /(XL(JFL+1)-XL(JFL)) 


C - - - CREATE  CAVITY  IN  ARRAYS 

C - (UPPER  SURFACE) 


I F (  TFU  .EQ.  1.0)  THEN 
JFU= JFU+ 1 
ELSE 

DO  90  I  =  KU , JFU+1 , - 1 
YU( 1+1 )=YU( I ) 

XU( l+1)=XU( I ) 

90  CONT I NUE 

C- . RESET  POINTERS 

KU=KU+ 1 
JFU= JFU+1 

c - INSERT  INTO  ARRAY 

YU(JFU)=YUFLAP 
XU(JFU)=AF 
END  IF 

C - (LOWER  SURFACE) 

IF(  TFL  .EQ.  1.0)  THEN 
JFL= JFL+ 1 
ELSE 

DO  100  I  =  KL , JFL+ 1,-1 
YL( 1+1 )=YL( I ) 

XL( 1+1 )=XL( I ) 

100  CONTINUE 

C - RESET  POINTERS 
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JFL= JFL+ 1 


C . . —  INSERT  INTO  ARRAY 

YL(JFL)=YLFLAP 
XL ( JFL ) =AF 
END  IF 


C - ROTATE  FLAP  ACCORDING  TO  DEFLECTION  ANGLE 

C  APPLY  COORDINATE  TRANSFORMATION 

C - UPPER  SURFACE 


XF=AF 

YF=YLFLAP+(YUFLAP-YLFLAP)/2 . 0 
FHYP=SQRT(XF*XF+YF*YF) 

FANG=ATAN2( YF ,  XF ) 

DO  170  I  =  JFU.KU 

CALL  FROTAT(XF ,  YF , FHYP , FANG, FDA . XU ( I ),YU( I ),XX,YY) 
XU( I )=XX 
YU ( I )=YY 
170  CONTINUE 

C - LOWER  SURFACE 

DO  180  I  =  JFL.KL 

CALL  FROTAT ( XF , YF , FHYP , FANG , FDA , XL ( I ) , YL ( I ) , XX , YY ) 
XL ( I )=XX 
YL( I )=YY 
180  CONTINUE 

I F (FDA  .GT.  0 . 0)THEN 


C - POSITIVE  FLAP  DEFLECTION 

C - FIND  HOW  MANY  POINTS  OVERLAP  ON  LOWER 

C  SURFACE 


I  FO=0 

DO  150  I  =  JFL.KL 

IF(XUJ)  .LT.  AF)  THEN 
I FO= I FO+ 1 


150 


END  IF 
CONT I NUE 


C - REMOVE  IFO  NUMBER  OF  OVERLAPPING  POINTS  FROM 

C  LOWER  SURFACE 

DO  160  I  =  JFL.KL-IFO 
XL ( I )=XL( I+IFO) 

YL( I )=YL( I+IFO) 

160  CONTINUE 

C - RESET  POINTERS 

KL=KL- I FO 

c - REMOVE  POINTS  FROM  UPPER  SURFACE  FOR 

C  CURVATURE  CONTINUITY 

IF  (TFU  .EQ.  1.0)  THEN 
NOPR= 1 


ELSE 

NOPR=2 
END  IF 

IF  (ABS(FDA)  .GT.  5*PI/180.)  THEN 
DO  183  I  =  JFU.KU-NOPR 
XU( I )=XU( l+NOPR) 

YU( I )=YU( l+NOPR) 

183  CONTINUE 

KU=KU-NOPR 


END  IF 
ELSE 

c - NEGATIVE  FLAP  DEFLECTION 

C . . FIND  HOW  MANY  POINTS  OVERLAP  ON  UPPER 

C  SURFACE 

I  F0=0 

DO  151  I  =  JFU.KU 


I F(XU( I )  .LT.  AF)  THEN 
I FO= I FO+ 1 
END  IF 


151  CONTINUE 

C - REMOVE  IFO  NUMBER  OF  OVERLAPPING  POINTS 

C  FROM  UPPER  SURFACE 

DO  161  I  =  JFU.KU- IFO 
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XU( I )=XU( I+IFO) 

YU( I )=YU( I+IFO) 

1 6 1  CONT I NUE 

C - RESET  POINTERS 

KU=KU- I FO 

C - REMOVE  POINTS  FROM  LOWER  SURFACE  FOR 

C  CURVATURE  CONTINUITY 


IF  (TFL  .EQ.  1.0)  THEN 
NOPR= 1 
ELSE 

NOPR=2 
END  IF 

IF  (ABS(FDA)  .GT.  5*PI/180.)  THEN 
DO  184  I  =  JFL , KL-NOPR 
XL ( I )=XL( l+NOPR) 

YL ( I )=YL( l+NOPR) 

184  CONTINUE 

KL=KL-NOPR 
END  IF 
END  IF 

6000  CONTINUE 

C - STORE  BACK  INTO  ORIGINAL  X,  Y  ARRAYS 

NBODY=0 

DO  200  I  =  KL,1,-1 
X(NBODY+1)=XL( I ) 

Y(NBODY+1)=YL( I ) 

NBODY=NBODY  + 1 
200  CONTINUE 

DO  190  I  =  2 , KU 

X(NBODY+1)=XU( I ) 

Y ( NBODY  + 1 ) = YU ( I ) 

NBODY=NBODY+1 
190  CONTINUE 

OPEN (2 ,F I LE= 'OUT1 .DAT ' ) 

DO  210  I  =  1, NBODY 
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WRITE(2,*)X( I ),Y( I  ) 

write(*,*)i,x(i),y(i) 

210  CONTINUE 
CLOSE(2) 

END 

C 

SUBPOUT  I NE  SROTAT ( XS , YS , SHYP , SANG , SDA , XP , YP , XX , YY ) 

XX=  XS  -  SHYP*COS ( SANG+SDA )  +  XP*COS(SDA)  -  YP*SIN(SDA) 
YY=  YS  -  SHYP*S I N(SANG+SDA)  +  XP*SIN(SDA)  +  YP*COS(SDA) 
RETURN 
END 
C 

SUBROUT  I NE  FROTAT ( XF , YF , FHYP , FANG , FDA , XP , YP , XX , YY ) 

XX=  XF  -  FHYP*COS ( -FANG+FDA )  +  XP*COS(FDA)  +  YP*SIN(FDA) 
YY=  YF  +  FHYP*S I N(-FANG+FDA)  -  XP*S1N(FDA)  +  YP*COS(FDA) 
RETURN 
END 
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Table  6 


Upper  surface  co-ordinates  Wortman  FX63-137 
airfoi 1 
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APPENDIX  E  :  Linearization  of  vorticity  and  stream 


function  equations 


The  vorticity  stream  function  equations  in  general 
coordinates  are. 

(stream  function) 
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6un 
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+  *>n  +  6ipn 

(42) 

where  n  and  n*i  are  successive  time  or  iteration  level 
substituting  equations  (42)  into  the  vorticity  equation  (41) 


4  J(un+  &>").+(*"  +  *pn  )„(«"+  6un  ).-(*>"  +  6Vn)_(on+  6u") 

^  n  t  s  n 

=  -k  (-5-  ("  ("% 

expanding  and  neglecting  terms  0(6) 


♦  *S*“i  ' 

“jjr, 


(43) 


Same  procedure  can  be  applied  to  the  stream  function  equation 
to  transform  it  into  the  6  (delta)  form,  equation  (40)  thus 
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becomes . 


it  should  be  noted  that  all  terms  appearing  as  changes 
(delta's  -  6)  between  successive  levels  are  unknown,  other 
exact  variables  are  assumed  known  from  the  previous  iteration 
level.  Thus  dropping  the  superscript  the  above  derived 
relations  take  the  final  form. 

J  It  +  W  V“e  ‘  uns*(  '  *t*V  -RsHrO'6"  W 

*  -  "5v  v«  *  "ss'+f"  “  >n>n  (45) 

(-5-  **£  •  Hr  %)t+  (-J-  6*n  -  -J-  {*£)n* 

*  -  (-5-  -  -J-  *n)  •  (-J-  ■  -J-  *e)  -  •**  (46) 
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APPENDIX  F 


Derivation  of  the  pressure  Poisson  equation 


The  non-dimens ional i zed  Navier  Stokes  equations  are 
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(47) 

(48) 

(49) 


Differentiating  equation  (48)  with  respect  to  x  and 
adding  to  equation  (49)  differentiated  with  respect  to  y 
results  in  the  following. 


|^(Eqn.  F . 2 )  +  §p(Eqn.  F.3)  s 


*  &} 


"15  {§72 (^]  '  afa?^)  +  "  §y  2 (4$ }  (50) 


the  above  expression  has  been  simplified  using  continuity 
(47).  Terms  or,  the  right  hand  side  (RHS)  can  be  expanded  as 
fol lows 
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d  r  du\ 

gj  [MgjJ  =  H  uxxx  ♦  Mx  uxx  +  dxx  ♦  „x  uxx 

5ydx(^5x)  “  ^  vxxy  +  ^x  vxy  +  ^xy  vx  +  ^y  vxx 

a2  r  3ul 

SySx^J  =  P  Uyyx  +  "x  uyy  +  **xy  uy  +  *V  uxy 

3y  (MSy)  ~  ^  vyyy  *  *V  vyy  +  Myy  vy  +  My  vyy 

grouping  and  rearranging  the  above  expanded  terms  results  in 
the  following  simplified  RHS  for  equation  (50) 
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substituting  the  above  expression  back  into  equation  (50) 
leads  to  the  following  Poisson  equation  for  pressure. 

_2n  foul 2  favl2  du  dv  .  i  /  _z„ 

7  P  =^J  +  N  +2  5P  33c  +  -Tfi  {  V  v 

+  Mx^Zu  +uxV2p  +  +  2  uy)}  (51) 


"h»r»  ^2=5y2  +  fl2  4 


using  the  transformat  ion  relations  given  in  Appendix  B  the 
Poisson  equation  for  pressure  in  general  coordinates  is 
obtained,  which  is: 
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The  above  relation  is  the  Poisson  equation  for  pressure 
in  general  coordinates  and  primitive  variables,  where  'o  ' 
represents  vorticity  and  ' u  * v  ’  the  velocity  components 
as  defined  in  the  physical  domain,  the  metrics  have  been 
defined  in  Appendix  B.  Note  that  all  terms  on  the  right  hand 
side  are  assumed  known  and  the  equations  is  solved  for  the 
only  unknown  variable  'P  ’pressure. 
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APPENDIX  Q  :  Derivation  of  force  and  moment  coefficients 

C  m  C  .  C 

l  a  «. 

c/4 

The  forces  acting  on  a  body  are  basically  the  body 
forces  (gravity)  the  shear  forces  (skin  friction)  and  the 
pressure  forces.  For  our  analysis  we  will  consider  only  the 
pressure  and  viscous  forces  (shear),  the  body  forces  are 
assumed  negilible  in  comparison  to  the  other  forces. 

Consider  the  dominant  forces  (pressure  and  shear)  acting 
on  the  differential  element  shown  in  the  following  figure 


Figure  34  Forces  acting  on  airfoil  element 


In  the  above  figure  x-y  are  the  body  axis  (fixed  to  the 
chord  of  the  airfoil)  and  X‘~Y’  are  the  wind  axis  (alligned 
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•nd  fixed  with  the  freestream).  The  aerodynamic  forces  are 

usually  measured  with  respect  to  the  wind  axis,  the  same  has 

been  followed  through  in  this  analysis. 

Fp  represents  the  force  due  to  pressure  acting  on  element 

i 

‘ i  '  and  is  defined  as 


Frr  pi  P«u.4sj 


(53) 


where  P.  is  non-dimensional  ized  pressure,  p  u  are  the  free 

1  CD  CD 

stream  properties  and  As .  is  the  element  area 

represents  the  force  due  to  shear  acting  on  element  ’ i  ' 


and  is  defined  as 


Ft  '  Cl  • 


p  u  As  . 


(54) 


where  is  the  skin  friction  coefficient  defined  a 


Using  principles  of  vector  analysis  the  forces  can  be 
decomposed  into  the  body  axis  system 


F  =  F.  cos  0.  +  F_  sin  0. 
x.  fi  i  i 

F  =  F.  sin  0.  +  F p  cos  0.  (56) 

7  i  v  j  r  i 

where  Fx  and  Fy  are  the  forces  along  the  body  axis  system, 
subscript  i  indicates  element  dependence  and  is  the  local 
slope,  given  by 


0 


(57) 
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ift  and  Drag  computation 


The  forces  in  the  body  axis  system  can  now  be 
transformed  into  lift  and  drag  using  the  following  relations 
reference  Figure  34. 


Li 
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F 
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F 

cos  a 

yi 

Xi 

where  L  .  and  D.  are  the  lift  and  drag  forces  acting  on 
element  i  respectively  and  a  is  the  airfoil  angle  of  attack. 
Sustituting  equations  (56)  into  the  relations  above  result 
in  the  following  simplified  form  for  equations  (58). 


L  .  =  F -  sin(O  .-a)  -  F_  cos (8  .-a) 

1  1  Fi  1 

D.  =  F-  cos(8.-a)  +  F  sin(8.-a)  (59) 

1  l  r i  1 

The  total  lift  and  drag  acting  on  the  body  is  simply  the 
summation  of  all  element  contributions,  thus 


M-l 

L  -  St, 


M-l 

D  =  I  Di 


(60) 


i-i  i-t 

where  n  is  the  number  of  nodes. 

Non-dimensional izing  the  above  relations  and  writing  in 
coefficient  form  equations  (60)  transform  to 


N-l 


C*  = 


— —  T  L  . 
o  u  c  i 

CD  oo  i-i 


°d  = 


p  u  c 
^00  00  1-1 


s-l 

i 


(61) 


where  and  now  represent  the  lift  and  drag  coefficients 
respectively  for  the  section.  Sustituting  equations  (59) 
into  the  above  relations  results  in  the  following  after 
simpl if i cat ion. 
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(62) 


Ct  =  £  |  Asi(c£  sinCe^-a)  -  Cp  cos^.-oojj 
Cd  =  E  |  As|-  [c^  cos(e.-a)  +  Cp  sin(0l.-a)j| 


Where  C_  =2  P  .  and  C-  =  -2  w/Re 

pi  1 


Moment  computation 

Traditionally  the  moment  acting  on  an  airfoil  is 

evaluated  at  the  quarter  chord  location,  and  is  therefore 

usually  known  as  the  "quarter  chord  moment".  To  evaluate 

this  quatity,  lets  consider  the  following  diagram,  where 

F  and  F  are  the  forces  acting  on  the  element  ’ i  '  given  by 
x  y 

relations  (56). 


Figure  35  Moment  acting  on  airfoil  element 


Then  the  moment  about  the  quarter  chord  (x  ,y  )  due  to  forces 

c  c 

acting  on  the  element  at  location  (xp,yp)  is  given  by: 
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(63) 


mc/4=  '  rTFt. 

where  is  the  radial  distance  and  F^  is  the  force  acting 

tangent  to  the  vector  r  .  (see  diagram  above). 

Decomposing  the  forces  F  and  F  into  the  radial  and 

a  y  • 

i  7 1 

tangential  components  results  in  the  following  expression  for 
the  tangential  force. 


F.  =  Fw  sin  0  +  F  cos  0 

t.  x.  -  y,- 


(64) 


substituting  relations  (56)  above  results  in  the  following 
after  simplification. 


F.  =  F-  cos (fl .  -  0  )  -  F  sin(0  .  +  0  ) 
ti  1  Pi  1 


(65) 


where  0^  is  given  by  equation  (57)  and  0^  is  defined  below 


0  =  tan 


1^) 


(66) 


substituting  equation  (65)  into  equation  (63)  and 
integrating  about  the  airfoil  results  in  the  following. 

*  "i’ {•-(  em)  -  F^sin(er  e.))j  (67) 


in  coefficient  form  the  above  relation  can  be  rewritten  as: 

(68) 


m 


=  .2  |ri  As^Cp  008(0^+  0^)  -  Sin(0.-  e.)]| 


where 


■( 


<*p  -  *  (yp  -  Ve)2)  2 


(69) 


109 


Appendix  H  :  Contour  plots;  Stream  Function,  Velocity  and 

yorfc1c1^y 


For  ell  cases  considered  in  this  study,  the  contour 
plots  have  been  plotted  for  the  stream  function,  velocity 
(x-component  of  the  total  velocity)  and  vorticity,  these  are 
given  as  Figures  36  to  51. 
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STREAM  FUNCTION 


STREAM  FUNCTION 


02H 


gure  40  Constant  stream  function,  velocity  and  vorticity 
contour  plots  for  HYortman  FX  63-137  airfoil, 
a  =  1 2°,  6i#=  0° ,  6tm=  0°,  Re  =  1x10l 
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VELOCITY 


STREAM  FUNCTION 


x/c 

Figure  51  Constant  strsam  function,  velocity  and  vortlclty 
contour  plots  for  Wortman  FX  63-137  airfoil, 

«  =  12°,  6  =  5°,  6.  =  20°,  Re  =  1x10s 

i*  t* 
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Appendix  I:  Variation  of  Pressure  and  Skin  Friction 
Coefficient  for  Wortman  FX  63-137  Airfoil 


For  all  cases  considered  in  this  study,  the  chordwise 
distribution  of  skin  friction  coefficient  (C^)  and  the 
pressure  coefficient  (Cp)  on  the  surface  of  the  Wortman  FX 
63-137  airfoil  has  been  plotted  vs.  the  chord  ratio  (x/c). 
same  are  attached  in  this  Appendix  as  Figures  52  through  67. 
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x/c 


Figure  53  Pressure  coefficient  (C  )  and  skin  friction 

P 

coefficient  (C^)  vs.  x/c 
a  =  0°,  0°,  6t>=  0°,  Re  =  1x10s 
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Figure  57  Pressure  coefficient  (cp)  *nd  ekin  friction 
coefficient  (C^)  vs.  x/c 
a  -  16°,  6=0°,  6  =  0°,  Re  =  IxlO6 

l*  t* 
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Figure  58  Preeeure  coefficient  (C  )  end  skin  friction 

P 

coefficient  (C^)  vs.  x/c 

o  =  -4°,  6=5°,  6  =  0°,  Re  =  1x10s 

i*  t« 
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Figure  59  Pressure  coefficient  (Cp)  *nd  skin  friction 
coefficient  (C^J  vs.  x/c 
a  =  0°,  6=5°,  8=0°,  Re  =  1x106 

!•  t« 
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Figure  60  Pressure  coefficient  (Cp)  *n<*  *kin  friction 
coefficient  (C,)  vs.  x/c 
°. 


...  /c  >  «nd  skin  friction 

re  61  Pressure  coefficient  (Cp) 

w  /r 
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